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Abstract.  This  paper  describes  an  algorithm  for  maintaining  an  approximating  triangu¬ 
lation  of  a  deforming  surface  in  R3.  The  surface  is  the  envelope  of  an  infinite  family  of 
spheres  defined  and  controlled  by  a  finite  collection  of  weighted  points.  The  triangulation 
adapts  dynamically  to  changing  shape,  curvature,  and  topology  of  the  surface. 


1.  Introduction 

This  paper  develops  a  fully  dynamic  algorithm  for  maintaining  a  triangulation  of  a  surface 
embedded  in  M3  that  changes  its  local  and  global  shape,  curvature,  and  topology  with 
time. 
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Motivation.  Deforming  surfaces  arise  in  moving  boundary  problems  of  physical  sim¬ 
ulation,  where  they  act  as  boundaries  of  spatial  domains  that  grow  and  shrink  with  time. 
An  example  is  the  boundary  between  the  solid  and  the  liquid  portions  of  metal  during 
solidification  [17].  Another  is  the  phase  boundary  in  a  solid  alloy  that  goes  through  the 
nucleation,  growth,  and  coarsening  stages  [1].  Moving  boundaries  also  arise  naturally 
in  mold  filling  processes,  both  for  metal  and  other  materials  [15].  Such  physical  pro¬ 
cesses  are  simulated  through  numerical  computations  facilitated  by  a  mesh  representing 
the  boundary  and/or  domain.  This  mesh  may  be  a  two-dimensional  triangulation  of  the 
surface,  or  a  three-dimensional  triangulation  of  space  on  one  or  both  sides  of  the  surface. 
The  numerical  methods  require  that  the  triangles  and  tetrahedra  used  in  the  triangulation 
be  well-shaped,  which  usually  means  they  have  small  aspect  ratio,  or,  equivalently,  they 
avoid  small  and  large  angles. 

We  are  also  interested  in  using  deforming  surfaces  in  the  modeling  of  molecules. 
Deformations  happen  naturally,  for  example  in  the  folding  process  of  proteins  [5],  Beyond 
natural  phenomena,  we  see  a  purpose  in  creating  artificial  deformations,  for  example  to 
interpolate  continuously  between  two  time-slices  in  a  molecular  dynamics  simulation, 
or  between  reconstructions  of  a  protein  from  two  different  crystallizations. 

Skin  Surfaces .  The  approach  to  deforming  surfaces  taken  in  this  paper  is  based  on  the 
technical  notion  of  skin  surfaces,  as  introduced  in  [8].  The  main  reason  for  this  choice 
is  the  existence  of  fast  combinatorial  algorithms  based  on  the  theory  of  alpha  shapes 
[9].  A  skin  surface  is  defined  by  a  finite  collection  of  spheres  in  R3.  We  can  think  of 
the  spheres  as  points  with  real  weights,  and  we  occasionally  prefer  this  interpretation 
to  avoid  confusion  with  the  various  other  types  of  spheres  that  arise  in  this  paper.  We 
derive  an  infinite  family  of  spheres  from  the  finite  collection  by  convex  combination 
and  shrinking.  The  skin  surface  is  the  envelope  of  this  family.  Even  though  the  family  is 
infinite,  the  surface  can  be  finitely  described  through  a  decomposition  into  a  collection 
of  quadratic  surface  patches.  Each  patch  is  the  portion  of  a  sphere  or  a  hyperboloid 
lying  inside  a  convex  polyhedron  obtained  by  shrinking  the  Minkowski  sum  of  corre¬ 
sponding  Delaunay  and  Voronoi  polyhedra.  In  each  case  the  sphere  or  hyperboloid  and 
the  containing  polyhedron  are  defined  by  k  <  4  weighted  points  (the  original  spheres). 
These  polyhedra  taken  together  form  a  finite  tiling  of  space,  which  we  refer  to  as  the 
mixed  complex.  The  correctness  of  the  essentially  combinatorial  surface  triangulation 
algorithm  relies  on  the  availability  of  exact  geometric  information,  possibly  in  symbolic 
form.  Most  important  in  this  context  is  the  maximum  curvature  at  a  given  surface  point, 
which  we  show  varies  continuously  over  the  surface  and  in  fact  can  be  extended  naturally 
to  a  continuous  function  throughout  space.  Equally  important  is  the  knowledge  about 
when,  where,  and  how  the  surface  changes  its  topological  type.  This  and  other  geomet¬ 
ric  information  is  readily  computable  from  the  dual  complex,  the  mixed  complex,  and 
the  decomposition  of  the  surface  defined  by  the  mixed  complex.  All  this  is  explained 
shortly. 

Triangulation.  For  computational  purposes  we  want  to  approximate  the  skin  surface 
by  a  two-dimensional  triangulation.  We  follow  the  convention  in  topology,  where  a 
triangulation  means  a  simplicial  complex  whose  underlying  space  is  homeomorphic  to 
the  surface.  The  triangulation  also  approximates  the  surface.  Specifically,  its  vertices 
lie  on  the  surface  and  their  spacing  depends  on  curvature.  The  algorithm  maintains  the 
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triangulation  through  local  restructuring  operations: 

•  it  moves  vertices  in  space  to  adapt  the  triangulation  to  changing  shape, 

•  it  adds  and  removes  vertices  to  adapt  the  local  density  to  the  local  maximum 
curvature, 

•  it  adjusts  connectivity  provided  by  edges  and  triangles  to  reflect  changing  topology. 

The  local  operations  are  automatic  and  follow  the  deformation  of  the  surface  dictated 
by  the  gradual  change  of  the  weighted  points  defining  it.  The  maximum  curvature  at 
each  surface  point  is  a  single  real  number,  so  our  adaptation  to  local  density  produces  an 
isotropic  triangulation.  We  gain  flexibility  by  permitting  the  triangles  to  deviate  somewhat 
from  the  equilateral  shape,  and  we  use  that  flexibility  for  obvious  geometric  reasons  but 
also  for  algorithmic  efficiency.  The  deviation  is  measured  as  circumradius  over  length 
of  the  shortest  edge,  and  the  algorithm  guarantees  that  this  ratio  never  exceeds  Q2/2. 
Here  Q  is  one  of  the  constants  on  which  the  algorithm  depends,  the  other  being  C,  which 
controls  how  well  the  triangulation  approximates  the  surface.  Q  controls  how  far  local 
density  can  deviate  from  strict  inverse  proportionality  to  local  maximum  curvature.  The 
two  constants  need  to  be  chosen  judiciously  in  order  to  guarantee  the  correctness  of  the 
algorithm. 

Outline.  The  technical  portion  of  this  paper  is  divided  into  three  parts  and  nine  sec¬ 
tions.  Part  I  provides  the  geometric  background.  It  consists  of  Section  2  describing  skin 
surfaces,  Section  3  showing  that  normal  direction  and  maximum  curvature  vary  slowly, 
and  Section  4  introducing  a  combinatorial  method  for  triangulating  the  surface.  Part  II 
explains  the  algorithm.  It  consists  of  Section  5  discussing  adaptation  to  changing  shape, 
Section  6  discussing  adaptation  to  changing  curvature,  and  Section  7  discussing  adapta¬ 
tion  to  changing  topology.  Part  HI  proves  the  algorithm  is  correct.  It  consists  of  Section  8 
analyzing  the  adaptation  to  curvature.  Section  9  detailing  the  various  operations  of  the 
algorithm,  and  Section  10  analyzing  the  adaptation  to  changing  topology.  Section  11 
concludes  the  paper  with  suggestions  for  future  work. 


Part  I.  Geometry 

The  three  sections  here  introduce  the  skin  surface,  analyze  its  tangent  and  curvature 
behavior,  and  show  that  with  a  dense  sampling  we  can  triangulate  the  surface  using  the 
restricted  Delaunay  triangulation. 


2.  Skin  Surfaces 

The  description  of  skin  surfaces  and  their  properties  offered  in  this  section  is  perhaps 
somewhat  terse.  The  reader  who  wishes  more  background  material  is  referred  to  [8]  for 
the  original  introduction  of  skin  surfaces,  to  [7]  and  [9]  for  a  description  of  alpha  shapes, 
and  to  [16]  for  a  textbook  in  geometry  that  talks  about  a  version  of  the  vector  space  of 
spheres  used  in  the  construction  of  skin  surfaces. 
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Sphere  Algebra.  Let  a  =  (a,  A)  be  the  sphere  with  center  a  e  R3  and  radius  A.  We 
require  A2  e  R.  For  A2  <  0  the  radius  is  imaginary  and  we  call  a  an  imaginary  sphere.  Its 

weighted  (square)  distance  function  7t a'  -+  R  is  defined  by  n&(x)  =  ||*  -a\\2-  A2; 

the  original  sphere  is  the  zero-set  of  this  function.  We  know  how  to  add  functions  and 
how  to  multiply  them  by  scalars.  If  we  apply  these  operations  to  the  7Xq  we  get  the  vector 
space  of  functions  of  the  form 

*(*)  =  Y  (11*  -  ptf  -  P), 

where  /5,  y  e  R  are  scalars  and  p  £  M3  is  a  point.  The  zero-set  of  n  is  the  sphere  with 
center  p  and  radius 

We  simplify  notation  by  applying  operations  directly  to  spheres.  In  particular,  we 
write  aa  +  fib  for  the  zero-set  of  an a  +  fin-b.  Using  this  notation,  we  can  define  what 
we  mean  by  the  affine  hull  and  by  the  convex  hull  of  a  finite  collection  of  spheres 
-4  =  [au  a2, . . . ,  a„],  namely 


aSA 


J2y‘a‘  =  1 

1=1  1  =  1 


conv.4  = 


T  Yt  €  aff  A  |  Yi  >  0  for  all  i  | 


As  an  exercise,  the  reader  may  want  to  verify  that  if  A  contains  only  two  spheres  and  they 
intersect  in  a  common  circle,  then  the  affine  hull  contains  exactly  all  spheres  passing 
through  this  circle.  The  convex  hull  contains  the  subset  whose  centers  lie  on  the  line 
segment  connecting  the  centers  of  the  two  given  spheres. 

Besides  adding  and  multiplying  with  a  scalar,  we  need  to  be  able  to  shrink  spheres. 
For  this  purpose  we  define  Vd  =  (■ a ,  A/V2),  which  is  the  zero-set  of  jra  +  A2/ 2. 
The  application  of  the  shrinking  operation  to  all  spheres  in  a  family  T  is  denoted  as 
s/p  =  {y/a  |  a  e  T).  The  skin  is  the  envelope  of  the  spheres  in  the  convex  hull  after 
shrinking, 

skin  A  —  env  VconvA. 

In  other  words,  the  skin  is  the  boundary  of  the  body,  denoted  body  A ,  which  is  the  union 
of  the  balls  bounded  by  spheres  in  Vconv  A 

Mixed  Cells.  The  mixed  cells  mentioned  in  the  Introduction  are  obtained  from  the 
corresponding  weighted  Voronoi  and  Delaunay  polyhedra.  For  a  given  finite  collection 
of  weighted  points  A ,  the  Voronoi  polyhedron  of  a  e  A  is  the  set  of  points  x  at  least  as 
close  to  a  as  to  any  other  weighted  point,  va  =  {*  €  R3  |  7ra(x)  <  for  all  b  e  A}. 
Two  Voronoi  polyhedra  meet  at  most  along  a  common  piece  of  their  boundary,  and  we 
define 

vx = n 

ieX 

for  every  subset  X  C  A.  It  is  convenient  to  assume  general  position,  in  which  case 
the  dimension  of  each  non-empty  v#  is  dim  vx  =  4  —  card  X.  In  particular,  vx  is  a 
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Fig.  1.  From  left  to  right:  a  typical  Voronoi  polyhedron,  a  Voronoi  polygon  times  a  Delaunay  edge,  a  Voronoi 
edge  times  a  Delaunay  triangle,  a  Delaunay  tetrahedron. 


polyhedron,  polygon,  edge,  vertex  if  X  has  cardinality  1,  2,  3,  4,  respectively.  Each 
non-empty  intersection  of  Voronoi  polyhedra  has  a  dual,  which  is  geometrically  realized 
as  the  convex  hull  of  the  (unweighted)  points  generating  the  polyhedra: 

8X  =  conv  {a  |  a  e  Xj. 

Assuming  general  position,  the  8#  are  simplices,  namely  vertices,  edges,  triangles,  tetra- 
hedra.  The  collection  of  these  simplices  is  referred  to  as  the  Delaunay  complex  of  A , 
although  usually  in  the  literature  this  term  is  reserved  for  the  case  of  unweighted  points. 

Note  that  vx  and  8X  have  complementary  dimensions:  dim  vx  +  dim  8X  =  3.  Fur¬ 
thermore,  they  lie  in  orthogonal  affine  subspaces  of  R3.  We  use  vector  operations  in  R3 
to  construct  the  mixed  cell  as  a  Minkowski  sum, 

fix  =  (Vx  +  <$Af)/2. 

The  dimension  of  fxx  is  always  3  =  dimv*  +  dim  8X.  Figure  1  shows  examples  of 
the  four  different  types  of  mixed  cells  corresponding  to  different  cardinalities  of  X.  The 
collection  of  mixed  cells  forms  a  face-to-face  tiling  of  R3,  which  we  call  the  mixed 
complex .  Figure  1 1  shows  the  mixed  complex  defined  by  four  points  in  the  plane. 

Skin  Patches.  Within  the  mixed  cell  fix,  the  skin  surface  is  completely  determined  by 
the  (at  most  four)  weighted  points  in  X  [8].  Specifically,  it  is  the  same  as  the  envelope 
of  the  affine  hull  after  shrinking  all  spheres,  that  is, 

skin *4  fl  fiX  =  envVaffA7  O  px. 

Let  k  =  card  A'  —  1.  Then  for  k  =  0  or  3  the  envelope  of  \f  dff  X  is  a  sphere,  and 
for  k  =  1  or  2  it  is  a  hyperboloid  of  revolution.  The  hyperboloids  have  asymptotic 
double-cones  with  right  opening  angles.  In  each  case  we  define  the  center  as  the  point 
zx  that  is  common  to  the  affine  subspaces  defined  by  vx  and  by  8X •  In  the  case  of  a 
hyperboloid  this  is  the  apex  of  the  asymptotic  double-cone,  and  in  the  case  of  a  sphere 
it  is  the  center.  It  may  or  may  not  belong  to  the  mixed  cell,  and  we  have  zx  £  fix  iff 
vx  O  Sx  is  non-empty. 

If  we  translate  the  center  to  the  origin  and,  in  the  case  of  a  hyperboloid  {k  =  1  or 
2),  rotate  so  that  the  axis  of  symmetry  is  along  the  x3-axis,  we  put  the  envelope  into 
standard  form.  If  R  is  the  minimum  distance  from  the  origin  to  the  envelope,  then  the 
equations  of  the  sphere  and  the  hyperboloids  are 

x\+x\+x  l  =  R2, 
x2+x2-x2  =  ±R2. 


(1) 

(2) 
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Fig.  2.  The  sphere,  the  one-sheeted  hyperboloid,  and  the  two-sheeted  hyperboloid. 


The  plus  sign  gives  the  one-sheeted  hyperboloid  and  the  minus  sign  gives  the  two-sheeted 
hyperboloid.  The  double-cone  arises  as  the  limiting  case  for  R  =  0.  The  three  surfaces 
are  illustrated  in  Fig.  2. 

Metamorphoses.  A  rather  simple  kind  of  deformation  of  the  skin  surface  is  generated 
by  increasing  the  weight  of  every  point  in  A  in  a  uniform  manner.  We  call  this  the  growth 
model  of  deformation  and  note  that  the  technical  results  in  this  paper  are  restricted  to 
this  model.  It  is  generated  by  changing  the  original  weights  A 2  of  the  weighted  points 
a  to  A2  -f  t  at  time  t.  It  is  easy  to  see  that  this  weight  change  preserves  the  Voronoi 
polyhedra  and  therefore  also  the  Delaunay  simplices  and  mixed  cells.  Even  though  the 
mixed  complex  remains  unaffected,  we  observe  all  generic  types  of  topological  changes 
or  metamorphoses  that  arise  in  general  deformations.  This  is  illustrated  in  Fig.  3.  As  indi¬ 
cated  in  Table  1 ,  there  are  four  types  depending  on  k  =  card  X  —  1 ,  where  px  is  the  mixed 
cell  containing  the  metamorphosis.  By  reversing  time  we  get  the  inverse  operations. 

We  can  also  reverse  the  orientation  of  the  skin  surface  by  finding  another  finite  collec¬ 
tion  of  weighted  points  that  has  the  same  skin  and  a  complementary  body.  Specifically,  we 
can  define  a  collection  of  spheres  B  =  A1  with  skin  A  =  skin  B  and  body  A  U  body  B  = 
M3.  Essentially,  B  contains  a  weighted  point  at  every  Voronoi  vertex  b  =  vx ,  and  the 
weight  is  chosen  so  that  B 2  =  || a  —b\\2  —  A 2  for  every  a  €  X  [8].  When  we  re¬ 
visit  the  metamorphoses  listed  in  Table  1  and  reinterpret  them  by  what  they  do  to  the 
body  of  B,  we  notice  a  symmetry  between  cases  k  and  3  —  k.  In  other  words,  there 
are  only  two  basic  types  of  metamorphoses.  The  first  type  is  geometrically  realized 


Fig.  3.  Three  snapshots  of  a  deforming  skin  triangulation  defined  by  continuously  growing  spheres.  From 
left  to  center,  note  in  the  front  two  k  =  1  metamorphoses,  each  adding  a  handle.  From  center  to  right,  note  at 
the  left  a  k  =  2  metamorphosis,  closing  a  tunnel. 
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Table  1.  The  four  types  of  generic  metamorphoses 
that  happen  during  growth/shrinking. 


k 

Type  of  metamorphosis/inverse 

0 

Creating/annihilating  a  component 

1 

Adding/removing  a  handle 

2 

Closing/opening  a  tunnel 

3 

Filling/starting  a  void 

by  a  sphere  appearing  or  disappearing.  The  limit  configuration  is  a  point,  and  in  the 
growth  model  this  is  the  center  of  the  appearing  or  disappearing  sphere.  The  second 
type  of  metamorphosis  is  geometrically  realized  by  a  two-sheeted  hyperboloid  flipping 
over  to  a  one-sheeted  hyperboloid,  or  vice  versa.  The  limit  configuration  is  a  double¬ 
cone,  and  in  the  growth  model  this  is  the  shared  asymptotic  double- cone  of  the  two 
hyperboloids. 

Time  of  Change.  An  interesting  question  is  when  exactly  the  metamorphoses  happen. 
We  answer  this  in  the  context  of  the  growth  model  by  introducing  certain  subcomplexes 
of  the  Delaunay  complex.  In  the  literature  these  subcomplexes  are  referred  to  as  dual 
or  alpha  complexes  [9],  but  we  use  different  notation  and  simply  denote  them  by  K  (t). 
Here  t  e  R  is  time  as  above.  The  growth  model  replaces  each  weight  A 2  by  A2  +  t  at 
time  t .  Restrict  each  Voronoi  polyhedron  to  within  the  generating  sphere  at  time  t,  giving 

Vfl(0  =  {*  e  Va  |  ||*  -  all2  <  A2  +  t}. 

The  complex  K(t)  consists  of  all  Delaunay  simplices  8x  for  which  the  restricted 
Voronoi  polyhedra  have  non-empty  intersection,  that  is,  C\azX  (0  7^  As  t  increases, 
K  —  K(t)  grows  into  a  progressively  larger  subcomplex  until  eventually  it  is  the  entire 
Delaunay  complex.  We  sort  the  simplices  in  the  order  they  enter  the  complex  K.  Even 
with  assumption  of  general  position  there  are  ties,  which  we  leave  unresolved,  by  allow¬ 
ing  more  than  one  simplex  at  a  given  position  in  the  ordering.  The  result  is  a  sequence 
of  collections  of  Delaunay  simplices  that  captures  the  evolution  of  the  complex.  Every 
prefix  of  the  sequence  is  itself  a  complex.  Because  of  this  property,  we  also  have  a  fast 
algorithm  for  deciding  how  and  when  the  homotopy  type  of  K  changes  [6]. 

The  underlying  space  of  K  (t)  and  the  body  bounded  by  the  skin  at  time  t  are  homotopy 
equivalent  [8].  It  follows  that  the  metamorphoses  for  the  two  structures  happen  at  exactly 
the  same  moments  in  time,  and  these  moments  can  be  computed  from  the  sequence  of 
simplices.  Assuming  general  position,  there  is  a  metamorphosis  for  every  position  in  the 
ordering  occupied  by  a  single  Delaunay  simplex  8 x ,  and  the  type  of  metamorphosis  is  the 
dimension  of  Sx  •  Whenever  there  are  two  or  more  simplices  tied  at  any  one  position,  their 
effects  on  the  homotopy  type  of  K  cancel  and  the  body  does  not  change  its  topological 
type. 

Sandwiching  Spheres.  We  close  this  section  by  stating  a  rather  special  and  important 
property  of  skin  surfaces,  from  [8],  which  is  heavily  exploited  in  this  paper.  We  mentioned 
already  that  the  skin  is  the  envelope  of  two  families  of  spheres,  one  inside  and  the  other 
outside  the  surface.  As  always  we  write  B  —  A 1. 
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Sandwich  Property.  For  every  point  x  on  the  skin  of  A,  there  are  unique  spheres 
Sx  e  V conv  A  and  Tx  e  Vconv  B  that  pass  through  x .  These  two  spheres  Sx  and  Tx  are 
externally  tangent,  and  have  equal  radius.  The  skin  surface  stays  outside  both  spheres, 
and  is  thus  tangent  to  them  at  x. 

We  refer  to  Sx  and  Tx  as  the  sandwiching  spheres  at  x  because  they  squeeze  the  surface 
fiat  in  a  neighborhood  of  x.  They  also  limit  the  normal  curvatures  at  x ,  and  we  will  see  in 
Section  3  that  they  in  fact  determine  the  maximum  curvature.  The  fact  that  Sx  and  Tx  are 
equally  large  follows  from  Lemma  7  in  [8].  In  a  nutshell,  the  reason  is  that  all  spheres 
S  6  conv  A  are  orthogonal  or  further  than  orthogonal  to  all  spheres  in  T  e  conv  B.  If 
we  shrink  S  and  T  each  to  y/2/2  of  the  original  size,  then  the  two  shrunken  spheres  are 
necessarily  disjoint,  unless  S  and  T  are  orthogonal  and  of  the  same  size,  in  which  case 
the  two  shrunken  spheres  touch  at  a  point.  Since  the  skin  is  the  common  envelope  of 
\! conv  A  and  \/conv  B,  the  two  spheres  passing  through  x  must  be  derived  from  equally 
large  spheres,  which  stay  equally  large  after  shrinking. 


3.  Continuity  of  Curvature 

This  section  proves  that  the  maximum  curvature  is  continuous  and  satisfies  a  Lipschitz 
condition.  We  use  this  to  control  local  density  in  the  triangulation.  This  section  also 
proves  a  one-sided  Lipschitz  condition  for  the  normal  direction. 

Maximum  Curvature.  Given  a  surface  F,  a  point  x  on  F,  and  a  tangent  vector  t*,  the 
normal  curvature  of  F  is  that  of  a  geodesic  passing  through  x  in  the  direction  tx.  The 
maximum  curvature  is  the  function  k:  F  R  that  maps  x  €  F  to  the  maximum  normal 
curvature  at  x.  For  a  hyperboloid  of  revolution,  the  minimum  curvature  is  measured 
within  planes  containing  the  symmetry  axis  (along  meridians),  and  the  maximum  cur¬ 
vature  is  measured  in  the  orthogonal  direction  (along  latitudes).  Explicit  expressions  for 
k  are  easy  to  compute  [12,  Chapter  14].  For  the  sphere  and  the  hyperboloids  in  standard 
form  (1)  and  (2),  the  maximum  curvatures  are 

*  =  1  /R,  (3) 

k  =  l/J±R2  +  2xl  (4) 

where  we  take  the  plus  sign  for  one-sheeted  hyperboloids  and  the  minus  sign  for  two- 
sheeted  hyperboloids.  By  plugging  ±R2  =  x2  +x%  — jef  into  (4)  we  see  that  the  maximum 
curvature  at  x  is  one  over  the  distance  of  x  from  the  origin.  This  implies  that  points  with 
constant  maximum  curvature  lie  on  spherical  shells  around  the  origin. 

Iso-curvature  Lemma.  Every  point  x  e  M3  belongs  to  exactly  one  hyperboloid  in 
standard  form,  andthe  maximum  curvature  of  that  hyperboloid  at  x  is  k{x)  =  1/||jc||. 

For  either  type  of  hyperboloid,  1//?  is  the  maximum  of  the  curvature  over  the  whole 
surface.  For  the  one-sheeted  hyperboloid,  R  is  also  the  radius  of  the  smallest  circle 
around  the  neck  of  the  hourglass.  For  the  two-sheeted  hyperboloid,  R  is  also  half  the 
smallest  distance  between  the  two  sheets. 
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Fig.  4.  The  circles  Sx  and  Tx  sandwich  the  hyperbola.  Depending  on  whether  we  revolve  the  hyperbola 
around  the  vertical  or  the  horizontal  axis,  we  get  a  one-sheeted  or  a  two-sheeted  hyperboloid. 


Curvature  Continuity.  To  prove  that  k  varies  continuously  over  the  skin  surface,  we 
consider  the  two  infinite  families  of  spheres  that  define  the  skin  as  their  common  envelope. 
For  a  finite  set  of  spheres  A,  let  S  =  V conv  A  and  T  =  V conv  A^-.  The  skin  of  A  is 
F  =  env  S  =  env  T.  The  family  S  defines  F  from  the  inside,  and  T  defines  it  from  the 
outside.  For  a  point  x  e  F,  there  are  unique  spheres  Sx  e  S  and  Tx  e  T  that  pass  through 
x.  We  make  essential  use  of  the  Sandwich  Property  stating  that  Sx  and  Tx  have  the  same 
size.  It  is  convenient  to  define  q(x)  =  1  /k(x)9  and  for  reasons  that  will  become  clear 
later  we  refer  to  q:  F  — >  E  as  the  length  scale. 

Curvature  Sandwich  Lemma.  For  every  point  x  €  F,  the  local  length  scale ,  g(x),  is 
the  common  radius  of  Sx  and  Tx . 


Proof.  If  x  belongs  to  a  sphere  patch,  then  that  patch  either  lies  on  Sx  or  on  Tx  and  q  O) 
is  obviously  the  radius.  Now  suppose  x  belongs  to  a  hyperboloid  patch.  The  hyperboloid 
is  obtained  by  revolving  a  hyperbola  around  one  of  its  two  symmetry  axes.  As  indicated 
in  Fig.  4,  the  hyperbola  is  the  common  envelope  of  two  families  of  circles,  one  centered 
along  each  of  the  two  symmetry  axes.  By  the  Sandwich  Property,  Sx  and  Tx  have  equal 
radii.  Because  x  is  halfway  between  the  centers  of  Sx  and  Tx,  that  radius  is  equal  to 
the  distance  of  x  from  the  origin.  By  the  Iso-curvature  Lemma,  this  distance  is  ||x|| 
=  q(x).  □ 


The  sandwiching  spheres,  and  their  common  radius,  vary  continuously  with  the  point 
x  €  F.  This  is  easy  to  see  for  points  in  the  interior  of  a  sphere  or  hyperboloid  patch, 
and  the  tangent  continuity  of  F  implies  the  same  for  points  on  the  boundary  common  to 
two  or  more  patches.  The  Curvature  Sandwich  Lemma  thus  implies  that  the  maximum 
curvature  varies  continuously  over  the  skin  surface  (except  at  centers,  where  it  blows 
up).  In  fact,  at  every  point  x  the  local  length  scale  q(x)  equals  the  distance  from  x  to  the 
center  zx  of  the  mixed  cell  [ix  that  contains  x. 
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Curvature  Variation.  We  strengthen  the  result  that  k(x)  is  continuous  by  showing  that 
it  varies  rather  slowly.  In  fact,  we  extend  its  reciprocal  £>(jc)  to  a  function  defined  on  all 
of  R3  and  show  that  £>(*)  has  Lipschitz  constant  one.  As  we  have  seen,  within  any  mixed 
cell  px,  Q  is  simply  the  distance  to  the  center  z  =  Zx-  By  the  definition  of  the  mixed 
complex,  this  is  a  continuous  function  on  R3.  Within  fix,  the  triangle  inequality  gives 
the  Lipschitz  bound, 

let*)  -  Q(y)\  =  III*  -zll  -  lly  -  zlll  <  II*  -  y\l 

By  applying  this  to  the  pieces  of  the  line  segment  from  x  to  y  contained  in  different 
mixed  cells,  we  obtain  the  result. 

Curvature  Variation  Lemma.  For  all  points  x,y  in  space  we  have  |e(*)  —  £(y)|  5 

ll*-yll. 

We  note  that  the  extension  of  q  to  a  function  R3  ->  R  describes  the  length  scale  of 
all  surfaces  in  the  family  defined  by  the  growth  model  of  deformation. 

Normal  Variation.  The  tangent  or  C1 -continuity  of  the  skin  surface  follows  from  the 
Sandwich  Property.  We  strengthen  this  result  by  proving  a  one-sided  Lipschitz  condition 
for  the  normal  vectors.  Specifically,  we  prove  an  upper  bound  that  relates  the  angle 
between  two  normal  vectors  at  points  x,  y  to  the  Euclidean  distance  between  them  and 
to  their  length  scales.  The  outward  unit  normal  vector  at  x  e  F  is  denoted  as  n*,  and 
the  angle  between  two  normals  is  Znxny  —  arccos(n^  •  ny).  In  proving  the  upper  bound, 
we  consider  again  the  one-parameter  family  of  skin  surfaces  generated  by  increasing 
square  radii  with  time.  For  points  x  =  (xi ,  x2,  *3)  on  a  sphere  in  standard  form  the  unit 
normals  are  nx  =  ±x/\\x\\,  and  for  points  x  on  a  hyperboloid  in  standard  form  they  are 
nx  =  ±(xi ,  x2 ,  — *3)/!!*  II  •  In  both  cases,  the  normals  are  the  same  along  a  line  passing 
through  the  origin,  and  they  vary  with  the  speed  of  the  angle  as  we  rotate  the  point  about 
the  origin.  The  formulas  imply  that  the  normals  of  points  x  and  y  in  two  adjacent  mixed 
cells  are  the  same  if  x  and  y  are  mirror  images  of  each  other  across  the  separating  plane. 
This  property  is  illustrated  in  Fig.  5. 
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Fig.  5.  Two  mixed  cells  with  dotted  circles  around  their  centers  and  parallel  normals  of  points  mirrored 
across  the  separating  edge.  The  illustration  shows  the  case  where  the  cells  have  k  —  0  and  1. 
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Normal  Variation  Lemma.  Let  x  and  y  be  points  on  F  with  distance  \\x  —  y||  <  £(*). 
The  angle  Znxny  between  the  surface  normals  at  x  and  y  is  at  most  arcsin  \\x  —  y\\/ 
q(x). 

Proof.  Consider  first  the  case  where  x  and  y  belong  to  the  same  mixed  cell,  and 
translate  the  coordinates  so  that  the  center  is  at  the  origin.  Given  x  and  the  distance 
\\x  —  y||  between  the  two  points,  the  angle  subtended  at  the  origin  is  a  maximum  if 
||jc H2  =  \\x  -  y ||2  +  ||y||2.  We  thus  have 

y  ^  .  \\x-y\\ 

Znxny  <  arcsin - — — , 

q(x) 

as  claimed. 

Consider  second  the  case  where  x  and  y  lie  in  different  mixed  cells.  The  directed 
line  segment  from  x  to  y  passes  through  i  >  1  planes  h\,hi, . . .  ,hi  separating  adjacent 
mixed  cells.  Let  pj  =  xy  n  hj  be  the  intersection  points  ordered  from  x  to  y.  We 
construct  a  polygonal  path  that  starts  at  x  and  whose  length  is  \\x  —  y\\.  It  is  obtained 
from  xy  by  reflecting  the  portion  after  pj  across  the  plane  hj,  for  j  =  i,  i  —  1, . . . ,  1 
in  this  order.  The  endpoint  y'  of  the  path  is  contained  inside  the  sphere  with  radius 
\\x  —  y  ||  around  x ,  which  implies  that  the  angle  between  x  and  y'  subtended  at  the  origin 
is  <p  <  arcsin(||x  —  y\\/g(x)).  Since  n-y  is  normal  to  the  sphere  or  hyperboloid  defined 
for  the  mixed  cell  of  x  that  passes  through  y\  (p  is  also  the  angle  between  nx  and  n^.  The 
claim  follows.  □ 

The  proof  of  the  Normal  Variation  Lemma  does  not  require  that  x  and  y  belong  to  the 
same  skin  surface.  The  claimed  inequality  holds  more  generally  for  any  points  x,  y  e  R3 
with  normals  defined  by  the  one-parameter  family  of  skin  surfaces  mentioned  above. 
Suppose  the  distance  between  x  and  y  is  \\x  —  y  ||  <  ££>(x).  Then  the  Normal  Variation 
Lemma  implies 

Zn^iVy  <  arcsine, 

which  is  the  form  used  most  often  in  this  paper. 


4.  Triangulation 

A  finite  set  V  c  F  is  an  e -sampling  if  for  every  point  x  e  F  there  is  a  vertex  a  e  V 
whose  distance  from  x  is  ||  jc  —  a  ||  <  eq(x).  The  goal  of  this  section  is  to  prove  that  the 
restricted  Delaunay  triangulation  defined  by  an  e-sampling  is  homeomorphic  to  the  skin 
surface,  provided  the  following  condition  holds: 

(I)  0  <  e  <  e0, 

where  eo  =  0.279 ...  is  a  root  of 

/  2e  \  2e 

/(e)  =  2cos  I  arcsin  - - 1-  arcsine  I  —  - — 

Note  that  /(e)  is  defined  for  —  1  <  e  <  and  that  it  is  non-negative  for  0  <  e  <  eo. 
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Restricted  Delaunay  Triangulation.  Let  V  be  a  finite  set  of  points  on  the  skin  surface. 
We  refer  to  these  points  as  vertices  and  we  denote  the  Voronoi  polyhedron  of  a  vertex 
a  €  V  in  K3  by  va.  The  corresponding  restricted  Voronoi  polygon  is  the  intersection  with 
the  skin  surface,  F  0  va,  which  is  non-empty  because  a  e  F  and  a  G  va.  The  restricted 
Delaunay  triangulation  is  the  nerve  of  the  collection  of  restricted  polygons: 


Dv  = 


conv U  \  U  C.  V,  Ffl  (~]va^0 

a<=U 


We  assume  general  position  and  in  particular  that  there  are  no  four  restricted  Voronoi 
polygons  with  non-empty  common  intersection.  It  follows  that  D  =  Dv  is  a  collection  of 
vertices,  edges,  and  triangles  but  contains  no  tetrahedra.  By  construction,  D  is  a  simplicial 
complex.  The  goal  of  this  section  is  to  prove  that,  for  e  satisfying  Condition  (I),  D  is 
a  triangulation  of  F.  Following  the  standard  topology  terminology  [2],  this  means  the 
underlying  space  of  D  is  homeomorphic  to  F.  As  shown  in  [10],  it  suffices  to  prove 
that  every  non-empty  common  intersection  of  restricted  Voronoi  polygons  is  a  closed 
topological  ball  of  the  appropriate  dimension,  namely  3  minus  the  number  of  polygons. 
If  this  is  the  case  we  say  D  has  the  closed  ball  property. 

We  formulate  this  property  in  terms  of  the  (unrestricted)  Voronoi  polyhedra.  By 
assumption  of  general  position,  the  intersection  of  k  +  1  =  2, 3, 4  Voronoi  polyhedra  is 
a  polygon,  edge,  vertex.  Depending  on  the  case,  the  intersection  with  the  skin  surface  is 
to  be 

case  k  =  0:  a  closed  disk, 
case  k  =  1:  empty  or  a  closed  interval, 
case  k  —  2:  empty  or  a  single  point, 
case  k  =  3:  empty. 

The  case  k  —  0  corresponds  to  a  single  Voronoi  polyhedron,  which  has  non-empty 
intersection  with  F  because  its  generating  point  lies  on  F.  We  establish  four  technical 
lemmas  in  preparation  of  proving  that  D  has  the  closed  ball  property. 

Distance  Claims .  If  two  surface  points  lie  in  the  same  Voronoi  polyhedron,  then  they 
cannot  be  far  from  each  other,  and  if  they  lie  on  a  line  that  is  almost  normal  to  the  surface, 
then  they  cannot  be  close  to  each  other.  We  quantify  the  first  claim  under  the  assumption 
that  V  is  an  ^-sampling. 


Short  Distance  Claim.  If  points  x  and  y  on  F  belong  to  a  common  Voronoi  polyhedron 
defined  by  a  vertex  in  an  e-sampling  V  c  F,  then  \\x  —y  ||  <  (2e/(l  —  e))g(x). 


Proof  Let  a  be  the  generating  point  of  the  common  Voronoi  polyhedron.  By  the  e- 
sampling  assumption  we  have  \\x  —  a ||  <  eg( x)  and  || y  —  a\\  <  eg(y).  Using  the 
triangle  inequality  we  get  \\x  —  y\\  <  e(g(x)  +  g(y)).  The  Curvature  Variation  Lemma 
now  implies 


Q(x)  >  Q(y)  -  \\x  y|| 

>  (1  -  e)g(y)  -  eg(x ), 
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and  hence  (1  +  s)q(x )  >  (1  —  s)g(y).  The  distance  between  x  and  y  is  therefore 


ll-K-yll  <  e 


(■♦£) 


e(x ) 


2s 

1  —  £ 


e(x), 


as  claimed.  □ 

We  get  a  better  bound  on  the  distance  if  one  of  the  points  generates  the  Voronoi 
polyhedron.  Assuming  x  —  a  we  get  ||a  —  y||  <  sg(y)  <  eg  (a)  +  s\\a  —  y||,  which 
implies 


lla-yll  <  7 —Q(a). 

1  —  £ 

We  need  this  version  of  the  Short  Distance  Claim  in  the  proof  of  the  Voronoi  Polyhedron 
Lemma  below. 

Next  we  quantify  the  second  claim,  which  is  independent  of  V. 

Long  Distance  Claim.  Suppose  a  line  meets  F  in  two  points  x  and  y  and  forms  an 
angle  smaller  than  £  with  the  surface  normal  at  x.  Then  \\x  —  y  ||  >  2 q(x)  cos£. 

Proof  By  the  Sandwich  Property,  there  are  two  spheres  of  radius  £(*)  that  both  pass 
through  x  and  locally  sandwich  the  surface.  The  line  meets  the  two  spheres  at  x  and 
at  points  at  distance  larger  than  2cos§£(jc)  on  both  sides.  The  skin  surface  contains 
no  points  inside  either  sandwiching  sphere,  which  implies  the  claimed  lower  bound  for 
||jc  —  y||.  □ 

We  play  the  Short  and  Long  Distance  Claims  against  each  other  and  thus  reach 
contradictions  proving  various  claims. 

Normal  Lemmas.  If  the  vertices  of  a  short  edge  or  a  triangle  with  small  circumcircle 
lie  on  the  skin  surface,  then  the  edge  or  triangle  lies  almost  flat.  We  quantify  both  claims. 
For  an  edge  ah  let  tab  =  (b  —  a)/\\h  -a  ||  be  the  unit  tangent  vector.  The  first  result  is 
an  immediate  corollary  to  the  Long  Distance  Claim: 

Edge  Normal  Lemma.  The  angle  between  an  edge  ab  and  the  surface  normal  at  its 
vertex  a  is  Z tabna  >  tt/2  —  arcsin(||a  —  b\\/2g(a)). 

A  common  use  of  the  Edge  Normal  Lemma  is  when  ab  belongs  to  the  restricted 
Delaunay  triangulation  of  an  ^-sampling.  Then  g(a)  >  (1  —  £)£(*),  where  x  is  a  point 
in  the  intersection  of  the  dual  Voronoi  polygon  with  the  skin  surface.  Hence  \\a  —  b\\  < 
2£q(x )  <  (2^/(1  -  e))£(a).  The  angle  between  ab  and  the  surface  normal  at  a  is  then 

7 r  .  e 

ZXabna  >  -  -  arcsin  - - . 

2  1  —  £ 
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Z  p(°)  ; 


Fig.  6.  The  dashed  sandwiching  spheres  meet  the  solid  sphere  around  a  in  two  parallel  dotted  circles.  Vertices 
b  and  c  are  placed  to  maximize  the  angle  between  the  triangle  normal  and  the  surface  normal  at  a . 


Next  we  consider  the  triangle  normal  lemma.  We  assume  the  angle  at  a  inside  the 
triangle  abc  is  no  smaller  than  the  angles  at  b  and  c.  Let  Rabc  be  the  radius  of  the 
circumcircle  and  let  nabc  be  the  outward  unit  normal  vector  of  abc . 

Triangle  Normal  Lemma.  If  a  is  a  vertex  of  the  triangle  abc  with  greatest  angle , 
then  the  angle  between  the  normal  of  abc  and  the  surface  normal  at  a  is  Z nabcna  < 
axcsm(2Rabc/Q(a)). 


Proof  Consider  the  two  spheres  of  radius  g(a)  that  locally  sandwich  the  surface  at  a,  as 
shown  in  Fig.  6.  The  face  angle  at  a  is  at  least  7r/3  and  the  length  of  the  edges  ab  and  ac 
is  at  most  2 Rabc  each.  To  compute  a  bound  on  the  angle  between  na  and  nabc  we  assume 
\\a  —  c||  <  || a  -  b\\  and  consider  the  sphere  with  radius  || a  —  b\\  around  a.  It  intersects 
the  sandwiching  spheres  in  two  parallel  circles.  Let  2X  be  the  distance  between  these 
two  circles  and  note  that  X/\\ a  —  b\\  =  || a  —  b\\/2g(a)  by  dropping  a  perpendicular 
from  z  to  the  midpoint  of  ab  and  using  similar  triangles.  Hence  2X  =  || a  —  b\\2/g(a). 
Since  the  angle  at  a  is  greater  than  or  equal  to  the  ones  at  b  and  c,  be  is  the  longest 
edge  of  abc.  The  angle  between  the  edge  be  and  the  planes  of  the  intersection  circles  is 
therefore  less  than 


arcsin 


II a  -  b\\2/Q(a) 


<  arcsin 


2  Rt 


abc 


\\b  —  c\\  """  g(a) 

This  is  an  upper  bound  for  the  angle  between  the  two  normal  vectors  at  a. 


□ 


Suppose  that  abc  belongs  to  the  restricted  Delaunay  triangulation  of  an  ^-sampling. 
Then  g(a)  >  ( 1  —  £)p(x),  where  x  is  a  point  of  the  intersection  between  the  dual  Voronoi 
edge  and  the  skin  surface.  Hence  Rabc  <  eg(x)  <  (6/(1  —  e))g(a).  The  angle  between 
the  two  normals  at  a  is  then 

/  ■  2e 
Znabcna  <  arcsin  - - . 

1  —  e 

Closed  Ball  Property.  We  are  now  ready  to  prove  the  closed  ball  property  for  the 
restricted  Delaunay  triangulation,  assuming  V  is  an  e-sampling  of  F  satisfying  Condi¬ 
tion  (I).  We  assume  general  position  and  consider  the  three  cases  in  turn:  first  Voronoi 
edges,  then  Voronoi  polygons,  and  finally  Voronoi  polyhedra. 
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Fig.  7.  A  Voronoi  polygon  intersecting  the  skin  in  a  circle  to  the  left  and  two  intervals  to  the  right. 


Voronoi  Edge  Lemma.  A  Voronoi  edge  of  V  intersects  the  skin  surface  in  at  most  one 
point. 

Proof.  Assume  there  is  a  Voronoi  edge  that  intersects  F  in  at  least  two  points,  *  and 
y .  Let  a  be  be  the  dual  triangle  in  the  restricted  Delaunay  triangulation.  The  Triangle 
Normal  Lemma  gives  an  upper  bound  for  the  angle  between  the  normal  of  abc  and  the 
surface  normal  at  a.  The  Normal  Variation  Lemma  gives  an  upper  bound  for  the  angle 
between  the  surface  normals  at  a  and  x.  Together  they  imply  an  upper  bound  for  the 
angle  f  between  nabc  and  nx  : 

£  <  Znabcna  +  Znanx 
.  2e 

<  arcsin - j-  aresm  s. 

1  -  £ 

The  angle  £  is  also  the  angle  between  the  Voronoi  edge  and  nx.  The  Long  Distance 
Claim  implies  \\x  —  y\\  >  2£>(jc)cos£,  which  by  Condition  (I)  contradicts  the  upper 
bound  || jc  —  y  ||  <  (2e/(l  —  £))£(»  implied  by  the  Short  Distance  Claim.  □ 

Voronoi  Polygon  Lemma.  The  intersection  of  a  Voronoi  polygon  of  V  with  the  skin 
surface  is  either  empty  or  a  closed  topological  interval. 

Proof.  Assume  there  is  a  Voronoi  polygon  whose  intersection  with  the  skin  surface 
contains  a  topological  circle  or  two  topological  intervals,  as  shown  in  Fig.  7.  Let  ab  be 
the  dual  edge  in  the  restricted  Delaunay  triangulation,  and  let  x  be  an  arbitrary  point  of 
the  intersection.  If  x  lies  on  a  circle,  then  let  L  be  the  line  in  the  plane  of  the  polygon 
that  intersects  the  circle  in  a  right  angle  at  x.  We  have  ZLnx  <  ZL'nx  for  any  line  V 
in  the  same  plane  and  passing  through  x.  Choose  LI  to  minimize  the  angle  with  na.  The 
Edge  Normal  Lemma  implies  an  upper  bound  for  the  angle  between  V  and  the  surface 
normal  at  a .  The  Normal  Variation  Lemma  implies  an  upper  bound  on  the  angle  between 
the  surface  normals  at  a  and  x.  Together  these  inequalities  imply 

£ 

ZLnx  <  arcsin - b  arcsin  s . 

1  -  £ 

This  angle  is  less  than  the  upper  bound  for  £  in  the  proof  of  the  Voronoi  Edge  Lemma, 
which  implies  a  contradiction  between  the  two  distance  claims. 
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In  the  case  of  two  intervals  let  L  be  a  line  connecting  x  to  the  closest  point  y  on  the 
other  interval.  If  y  lies  in  the  interior,  then  L  intersects  the  interval  in  a  right  angle  at 
y .  In  this  case  we  get  a  contradiction  with  the  same  argument  as  above  only  with  x  and 
y  interchanged.  Otherwise,  y  is  an  endpoint  of  the  interval  and  lies  on  a  Voronoi  edge. 
The  angle  between  L  and  is  less  than  that  between  the  Voronoi  edge  and  ny.  We  thus 
get  a  contradiction  with  the  same  argument  as  used  in  the  proof  of  the  Voronoi  Edge 
Lemma.  □ 

Voronoi  Polyhedron  Lemma.  The  intersection  of  a  Voronoi  polyhedra  ofV  with  the 
skin  surface  is  a  closed  topological  disk. 

Proof  Assume  there  is  a  Voronoi  polyhedron  whose  intersection  with  the  skin  surface 
contains  a  closed  2-manifold  (without  boundary),  a  2-manifold  with  boundary  other  than 
a  disk,  or  two  disks.  In  the  first  case  we  let  L  be  a  line  that  intersects  the  2-manifold  in 
two  points,  x  and  y,  and  forms  a  right  angle  at  x.  We  get  a  contradiction  between  the 
two  distance  claims  as  before. 

For  the  rest  of  the  proof,  let  a  be  the  generating  vertex  of  the  Voronoi  polyhedron  and 
assume  the  intersection  between  this  polyhedron  and  the  skin  surface  is  a  2-manifold 
with  boundary,  F'.  This  2-manifold  with  boundary  can  be  different  from  a  disk  either 
because  it  is  non-orientable,  it  contains  a  handle,  or  it  has  at  least  two  boundary  circles. 
The  non-orientability  of  F'  contradicts  the  orientability  of  F.  If  Ff  has  a  handle  but  only 
one  boundary  circle,  then  homology  theory  gives  us  a  pair  of  simple  closed  curves  in  Fr 
that  intersect  each  other  transversely  exactly  once.  Along  either  one  of  these  curves,  there 
is  a  point  such  that  the  line  normal  to  Ff  that  passes  through  that  point  meets  the  other 
curve,  and  hence  Ff  again.  This  gives  a  contradiction  to  the  two  distance  claims.  A  more 
elaborate  argument  is  needed  for  the  case  where  there  are  two  or  more  boundary  circles. 
Then  either  F'  is  connected,  and  in  the  simplest  case  is  an  annulus,  or  it  is  disconnected, 
and  in  the  simplest  case  consists  of  two  disks. 

By  the  remark  after  the  Short  Distance  Claim,  the  distance  between  a  and  a  point 
y  €  Ff  is  || a  -  y||  <  ( e/(\  —  e))Q{a).  Let  L  be  the  normal  line  at  a  and  note  that  it 
contains  the  line  segment  of  length  2q (a)  that  connects  the  centers  of  the  two  spheres 
sandwiching  the  surface  at  a .  This  line  segment  is  contained  in  the  Voronoi  polyhedron, 
which  implies  that  the  polyhedron  is  fairly  tall  and  slim.  Consider  a  plane  that  contains  L 
and  intersects  at  least  two  boundary  circles  of  Ff .  Such  a  plane  exists  for  else  we  can  find 
a  plane  through  L  that  intersects  no  boundary  circle  at  all.  However,  then  L  meets  F'  in  at 
least  two  points,  and  we  again  get  a  contradiction  to  the  two  distance  lemmas.  The  plane 
that  meets  two  boundary  circles  intersects  the  Voronoi  polyhedron  in  a  convex  polygon 
and  Ff  in  at  least  two  connected  curves.  One  of  the  curves  contains  a.  We  may  assume 
that  the  second  curve  lies  on  one  side  of  L.  Let  V  be  the  line  passing  through  its  two 
endpoints,  which  both  lie  on  the  boundary  of  the  convex  polygon.  The  line  V  intersects 
the  sphere  with  radius  (s/(l  —  £))£(a)  around  a  and  it  does  not  intersect  the  line  segment 
connecting  the  centers  of  the  two  sandwiching  spheres.  The  angle  between  V  and  the 
surface  normal  at  a  is  therefore  ZnflL'  <  arcsin(f/(l  —  s)).  By  the  intermediate  value 
theorem  there  is  a  point  y  on  the  second  curve  whose  curve  normal  n'^  is  also  normal  to 
Z/.  Hence  Znfln^  >  tt/2  —  arcsin(£/(l  —  £)).  Since  the  surface  normal  at  y  is  also  normal 
to  the  tangent  line  parallel  to  L\  its  angle  with  na  is  at  least  this  large.  From  the  Normal 
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Variation  Lemma  we  get  Z.nany  <  arcsin(£/(l  -  e)).  Putting  both  inequalities  together 
we  get  7t/2  <  2arcsin(£/(l  —  £)).  This  is  equivalent  to  e  >  \/2  —  1  =  0.414 . . .  and 
contradicts  Condition  (I).  This  completes  the  proof  of  the  Voronoi  Polyhedron  Lemma 
for  the  final  case  where  F'  has  at  least  two  boundary  circles.  □ 

Summary.  The  three  Voronoi  lemmas  establish  that  for  e  satisfying  Condition  (I),  the 
restricted  Voronoi  diagram  of  an  ^-sampling  V  has  the  closed  ball  property.  The  result 
of  [10]  implies  that  the  underlying  space  of  the  restricted  Delaunay  triangulation  is 
homeomorphic  to  the  skin  surface. 

General  Homeomorphism  Theorem.  The  restricted  Delaunay  triangulation  of  an  e- 
sampling  triangulates  the  skin  surface,  for  e  satisfying  Condition  (I). 

For  the  purpose  of  changing  the  topology  of  the  skin  surface  we  rely  on  point  dis¬ 
tributions  that  locally  violate  the  e-sampling  condition.  We  give  a  separate  proof  of  the 
closed  ball  property  in  Section  10  and  thus  obtain  a  Special  Homeomorphism  Theorem 
for  such  distributions. 


Part  n.  Algorithm 

The  algorithm  maintains  the  triangulation  of  a  deforming  skin  surface  dynamically  by 
adapting  geometric  position  to  shape,  density  to  curvature,  and  connectivity  to  topology. 
It  can  be  used  to  construct  a  triangulation  by  starting  with  the  empty  triangulation  and 
growing  components  from  nothing. 


5.  Shape  Adaptation 

This  section  describes  the  overall  algorithm  and  presents  the  details  for  adapting  the 
triangulation  to  the  changing  shape  of  the  surface.  We  restrict  the  deformation  to  the 
growth  model,  where  the  weight  A2  of  every  sphere  is  changed  to  A2  +  t  at  time  t.  Let 
to  <  h  be  moments  in  time  and  let  D0,  A  be  the  corresponding  restricted  Delaunay 
triangulations.  The  algorithm  updates  Do  locally  and  changes  it  to  D\. 

Moving  Vertices .  The  intuition  for  moving  vertices  is  taken  from  Morse  theory,  which 
considers  structures  that  arise  in  sweeping  out  a  smooth  manifold  [14].  The  skin  surface 
is  the  cross  section  at  a  moment  in  time  during  the  sweep,  and  the  manifold  is  the  stack 
of  cross  sections  in  the  time  direction.  In  other  words,  the  manifold  is  the  graph  of 
M:  R3  — ►  R  that  maps  a  point  x  to  the  time  t  at  which  x  belongs  to  the  surface  F(t). 
Hence  F(t)  =  M~l  (t).  A  metamorphosis  of  F  corresponds  to  a  critical  point  of  M.  For 
cross  sections  in  a  time  interval  [to,  h]  that  is  free  of  critical  points,  we  can  construct  a 
one-parameter  family  of  diffeomorphisms  from  the  integral  lines  of  the  gradient  vector 
field  grad  Af(jc).  These  diffeomorphisms  <p,*:  F(t0)  -►  F(tt),  with  U  e  [t0,  h],  can  be 
composed  to  connect  diffeomorphically  any  two  cross  sections  in  the  time  interval, 

<Pij  =  <Pj  O  <pp:  F(ti )  ->  F(tj). 
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Fig.  8.  Dotted  integral  lines  of  a  solid  growing  circle  and  a  solid  growing  hyperbola. 


The  step  from  time  f,*  to  time  tj  thus  amounts  to  moving  each  vertex  a  e  D,  along  its 
integral  line  to  a'  —  <Pij(a)  e  Dj.  In  the  growth  model  the  integral  lines  are  pieces  of 
straight  lines  and  hyperbolas,  as  illustrated  in  Fig.  8.  To  see  this,  note  that  {lx\ ,  —2x3)  are 
the  normal  vectors  ofthe  family  of  hyperbolas  =  ±R2,  and  that  (2jc3,  2x\)  are  the 

normal  vectors  of  the  family  2*i*3  =  ±/?2  obtained  by  rotating  the  first  family  through 
an  angle  of  7r/4.  The  three-dimensional  picture  is  obtained  by  revolving  the  hyperbolas 
in  Fig.  8  about  the  x3-axis.  The  first  family  of  hyperbolas  turns  into  the  one-parameter 
family  of  hyperboloids  described  by  (2).  The  second  family  turns  into  a  two-parameter 
family  of  hyperbolas  each  orthogonal  to  each  of  these  hyperboloids. 

For  deformations  more  general  than  the  ones  in  the  growth  model,  we  may  not  be 
able  to  determine  the  integral  lines  explicitly.  Fortunately,  moving  vertices  along  integral 
lines  is  convenient  but  not  necessary  for  the  algorithm,  and  an  approximation  of  that 
movement  will  in  general  suffice.  For  small  time  steps,  the  triangulation  changes  only  a 
small  amount  and  can  be  maintained  with  the  methods  described  in  this  and  the  following 
two  sections. 

Parametrization.  It  is  convenient  to  parametrize  the  integral  lines  by  time  so  that  points 
can  be  moved  by  evaluation.  Each  integral  line  is  decomposed  by  the  mixed  complex 
into  pieces  of  lines  and  hyperbolas.  We  first  consider  the  case  of  a  line  inside  a  mixed  cell 
constructed  from  a  Delaunay  vertex  and  its  dual  Voronoi  polyhedron.  After  translating  the 
center  to  the  origin,  the  mixed  cell  is  swept  out  by  a  sphere  in  standard  form  = 

s,  for  s  >  0.  We  thus  get  integral  lines  that  start  at  the  origin  and  go  to  infinity,  and  we 
clip  each  such  half-line  to  within  the  mixed  cell.  If  the  origin  lies  inside  the  mixed  cell, 
then  it  is  the  source  of  an  entire  sphere  of  integral  lines.  We  follow  the  usual  convention 
and  parametrize  that  sphere  by  longitude  and  latitude,  0  e  [0,  27t)  and  <p  €  [ — tt,  it].  For 
each  pair  of  angles  we  have  a  half-line  y0y.  R+  U  {0}  ->  R3  defined  by 

(cos  0  cos  (p 
sin  0  cos  <p  ^/s  I  . 
sin  <p  *Js  ) 

The  case  of  a  mixed  cell  constructed  from  a  Delaunay  tetrahedron  and  its  dual  Voronoi 
vertex  is  symmetric,  with  the  integral  lines  ending  rather  than  starting  at  the  origin.  If  the 
origin  lies  inside  the  mixed  cell,  then  it  is  the  sink  of  an  entire  sphere  of  integral  lines. 

We  next  consider  the  case  of  a  mixed  cell  constructed  from  a  Delaunay  edge  and  its 
dual  Voronoi  polygon.  We  assume  the  hyperboloid  sweeping  out  the  mixed  cell  is  in 
standard  form  x\  +x\  —  x\  =  s,  for  s  e  R.  The  integral  lines  are  hyperbolas  that  fall  and 
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rise  along  the  jt3-axis  and  turn  a  total  of  90°  before  reaching  the  *i*2-plane,  which  they 
approach  as  they  disappear  to  infinity.  We  parametrize  the  family  with  the  longitudinal 
angle,  0  e  [0,  27r),  and  the  minimum  distance  to  the  origin,  R  >  0.  For  each  pair  of 
parameters  we  get  a  hyperbola  Yq,r'-  M  — >  M3  defined  by 

(±  cos  0  */u\ 

±sin0v^  J  , 


with  u  —  (s  4-aA2  4-  R4)  /2.  To  check  the  correctness  of  the  parametrization  note  that  the 
points  yotR(s)  satisfy  the  equation  of  the  hyperboloid  and  the  equations  of  the  orthogonal 
hyperbola.  The  case  of  a  mixed  cell  constructed  from  a  Delaunay  triangle  and  its  dual 
Voronoi  edge  is  symmetric,  with  the  integral  lines  moving  in  above  and  below  the  X\X2- 
plane  and  turning  a  total  of  90°  before  reaching  the  Jt3-axis,  which  they  approach  as  they 
go  to  infinity. 

In  either  case  we  obtain  a  parametrization  in  time  by  setting  s  =  R2  +  t.  Note  that 
in  all  four  cases  of  integral  lines,  the  speed  of  the  parametrization  depends  only  on  the 
distance  to  the  center  of  the  mixed  cell,  \\dy/ds\\  =  1/(2  ||y(s)||).  This  is  consistent 
with  the  length  of  the  gradient  of  M{x)  =  ±x\  ±x\±x\  being  independent  of  the 
choice  of  signs,  ||grad  Af(x)||  =  2||*||. 

Algorithmic  Time-Warp.  Vertices  move  continuously  along  their  integral  lines,  but 
updating  them  continuously  is  computationally  infeasible.  The  common  escape  from 
this  dilemma  is  the  time-slicing  method,  which  takes  discrete  time  steps  and  advances 
all  vertices  from  time  t0  to  time  tx  without  intermediate  stop.  There  are  drawbacks  to  time¬ 
slicing  related  to  the  difficulty  of  choosing  the  right  step  size.  We  follow  an  alternative 
approach  and  take  different  time  steps  at  different  locations.  This  is  done  by  prioritizing 
the  four  types  of  operations  that  occur  at  discrete  moments  in  time,  which  are  edge  flips, 
edge  contractions,  vertex  insertions,  and  metamorphoses.  Edge  flips  are  described  below. 
Edge  contractions  and  vertex  insertions  arise  in  curvature  adaptation,  and  are  described  in 
the  next  section,  while  metamorphoses  are  the  operations  that  allow  topology  adaptation. 

Coordinate  updates  are  done  lazily,  moving  a  vertex  when  and  only  when  it  is  used 
by  one  of  the  other  four  operations.  This  results  in  a  time-warped  surface  with  different 
pieces  reflecting  the  state  at  different  times.  To  bring  the  entire  surface  to  the  present 
time,  we  simply  update  all  the  vertex  coordinates,  and  by  assumed  correctness  of  the 
prioritization  this  requires  no  other  changes  in  the  triangulation. 

At  any  moment  in  time  t ,  we  consider  the  collection  of  possible  next  operations.  Let 
ti  >  t  be  the  time  at  which  such  an  operation  xt  would  happen  if  the  vertices  moved 
along  integral  lines  and  no  other  operations  preceded  .  We  store  the  r ,  in  a  priority 
queue  ordered  by  time.  The  overall  algorithm  is  a  simple  infinite  loop: 

loop  r i  =  NextOp;  D  =  Apply(t,)  forever. 

Function  Apply  changes  D  according  to  r,- ,  and  simultaneously  updates  the  priority  queue 
by  inserting  new  operations  made  possible  by  the  changes  caused  by  t*  .  The  changes  may 
make  some  of  the  operations  in  the  priority  queue  inapplicable.  For  example,  the  edge  of 
an  edge  flip  may  disappear  from  D.  Instead  of  deleting  these  operations  immediately,  we 
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Fig.  9.  Head-on  view  of  the  edge  be.  The  dotted  line  represents  the  skin  and  the  shading  indicates  the  side 
of  the  body.  To  the  left  be  is  convex  and  to  the  right  it  is  concave. 


use  a  lazy  strategy  that  checks  an  operation  when  it  reaches  the  top  of  the  priority  queue. 

Operation  NextOP: 

repeat  r  =  EXTRACTMlN  until  ISOk(t); 
return  z. 

Determining  when  exactly  an  operation  r,  matures  in  the  future  is  computationally  fairly 
expensive,  and  so  is  the  correct  ordering  of  operations  in  time.  We  plan  to  discuss 
approximate  ordering  methods  that  alleviate  the  cost  in  a  later  paper. 

Edge  Flipping.  Let  be  be  an  edge  of  the  restricted  Delaunay  triangulation  D  at  time 
t.  It  is  shared  by  two  triangles,  abc  and  bed.  By  the  Voronoi  Edge  Lemma,  the  Voronoi 
edges  dual  to  abc  and  bed  meet  the  skin  surface  in  a  point  each.  Let  La  and  Ld  be 
the  lines  that  contain  the  two  Voronoi  edges  and  orient  them  from  where  they  meet  F 
towards  the  point  x  =  La  O  Lj  where  they  cross.  The  point  x  may  or  may  not  be  a 
Voronoi  vertex.  Call  be  convex  or  concave  depending  on  whether  the  dihedral  angle 
between  abc  and  bed  measured  on  the  side  of  the  body  is  less  than  or  greater  than  tt. 
As  illustrated  in  Fig.  9,  in  the  convex  case  the  two  lines  pass  from  outside  to  inside  the 
body,  and  in  the  concave  case  it  is  the  other  way  round.  Flipping  the  edge  be  means 
replacing  it  by  the  other  diagonal  of  the  quadrangle  be  defines.  The  operation  can  be 
performed  unless  ad  is  already  an  edge  in  the  triangulation,  in  which  case  either  b  or  c 
belongs  to  only  three  edges.  The  flip  would  then  decrease  that  number  to  two  edges  and 
contradict  the  closed  ball  property  of  the  restricted  Delaunay  triangulation.  The  three 
Voronoi  lemmas  thus  imply  that  the  flip  of  be  would  not  be  attempted  if  ad  is  already 
in  the  triangulation. 


void  EdgeFlip(£c): 
assert  ad  g  D\ 

substitute  cad ,  ad ,  a db  for  abc ,  be ,  bed. 

The  edge  flip  operation  is  a  response  to  the  event  that  the  Voronoi  edges  dual  to  abc  and 
bed  stop  meeting  the  skin  surface.  This  happens  when  x  passes  through  F,  and  in  this 
case  x  is  necessarily  a  vertex  of  the  Voronoi  diagram.  If  be  is  convex,  then  x  passes  from 
inside  to  outside  the  body  and  ad  is  concave.  Symmetrically,  if  be  is  concave,  then  x 
passes  from  outside  to  inside  the  body  and  ad  is  convex.  The  time  f,  when  the  flip  happens 
depends  on  the  points  a,  b,  c,  d  and  the  surface  F,  all  of  which  move  continuously  with 
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time.  In  other  words,  f,-  is  a  root  of  a  continuous  function  in  t .  It  is  used  as  the  priority  of 
the  edge  flip  r *  stored  in  the  event  scheduling  priority  queue. 


6.  Curvature  Adaptation 

This  section  focuses  on  the  density  adaptation  algorithm,  implemented  through  edge 
contraction  and  vertex  insertion.  The  method  is  straightforward,  but  we  need  some 
geometric  analysis  to  convince  ourselves  that  it  is  correct. 

Invariants.  The  goal  of  the  algorithm  is  to  triangulate  locally  with  edges  and  triangles 
of  size  roughly  proportional  to  the  length  scale,  which  Section  3  defined  as  one  over  the 
maximum  curvature,  £(x)  =  1  /k(x).  The  size  of  an  edge  ab  is  defined  to  be  half  its 
length,  Rab  =  \\a  —  b\\/2.  The  size  of  a  triangle  abc  is  the  radius  Rabc  of  the  circumcircle, 
which  length  scale  they  should  follow  exactly.  For  edges  we  worry  about  them  getting 
too  short,  so  we  compare  size  with  the  maximum  length  scale,  and  for  triangles  we  worry 
about  them  growing  too  large,  so  we  compare  size  with  the  minimum  length  scale: 

Qab  =  ma x{g(a),Q(b)}, 
gabc  =  min{e(a),  gib),  q(c)}. 

The  algorithm  is  formulated  using  two  positive  constants,  C  and  Q .  Roughly,  C  controls 
how  closely  the  triangulation  approximates  the  skin  surface,  and  Q  controls  the  quality 
of  the  triangles.  The  following  two  inequalities  are  maintained  as  invariants,  which  we 
refer  to  as  the  Lower  Size  Bound  and  the  Upper  Size  Bound: 

[L]  Rab /Qab  >  C/Q  for  every  edge  ab  e  D. 

[U]  Rabc/Qabc  <  C  Q  for  every  triangle  abc  e  D. 

It  is  not  necessary  to  check  for  long  edges  and  small  triangles  explicitly.  This  is  because  an 
edge  of  size  Rab  >CQgab  belongs  to  two  triangles  that  both  violate  [U],  Symmetrically, 
a  triangle  of  size  Rabc  <  ( C/Q)Qabc  has  three  edges  that  violate  [L].  Appropriate  values 
of  C,  Q  will  be  determined  in  the  analysis  of  the  algorithm  but  we  can  already  anticipate 
C  =  0.08,  Q  =  1.65  as  a  feasible  assignment. 

Minimum  Angle.  The  smallest  angle  is  a  measure  of  triangle  quality.  It  achieves  its 
maximum,  n/ 3,  for  the  equilateral  triangle.  Triangles  that  satisfy  both  Size  Bounds 
cannot  have  arbitrarily  small  angles. 

Minimum  Angle  Lemma.  A  triangle  that  satisfies  [L]  and  [U]  has  minimum  angle 
larger  than  arcsin(l/  Q2). 

Proof.  Let  abc  be  a  triangle,  with  be  its  shortest  edge  and  R  its  circumradius.  We  get 
Qabcf  Qbc  <  1  by  definition  of  length  scale.  Using  [L]  and  [U]  we  get 

R  CQeabc  <  Q? 

\\b  -c II  <  2(C/ Q)Qbc  -  2  • 

The  minimum  angle  is  Zbac ,  and  \\b  —  c\\  =  2R  sin  Zbac.  Hence  Zbac  =  arcsin 
(|| b  —  c\\/2R)  >  arcsin(l/<22)»  as  claimed.  □ 


546 


H.-L.  Cheng,  T.  K.  Dey,  H.  Edelsbrunner,  and  J.  Sullivan 


The  Minimum  Angle  Lemma  suggests  that  we  choose  Q  as  small  as  possible,  contin¬ 
gent  upon  satisfying  all  constraints  needed  to  prove  the  algorithm  correct.  For  Q  =  1.65 
the  minimum  angle  is  larger  than  21.54 . . .°,  and  the  maximum  angle  is  smaller  than 
180°  -2*21.54°  =  136.90...°. 

Enforcement.  The  algorithm  enforces  the  two  invariants  by  contracting  short  edges  and 
inserting  vertices  near  the  barycenters  of  large  triangles.  Let  abc  be  a  triangle  that  gets 
too  large,  that  is,  Rabc  =  C Qgabc  at  time  The  time  r,  depends  on  the  points  a,  b,  c  and 
their  length  scales  g(a),  g(b ),  p(c),  which  all  change  continuously  with  time.  To  remedy 
the  violation  of  the  Upper  Size  Bound,  we  add  the  restricted  Voronoi  vertex  x  dual  to 
abc  as  a  new  vertex  to  the  tri angulation.  A  vertex  insertion  may  cause  new  violations  of 
the  Upper  Size  Bound  and  thus  trigger  additional  vertex  insertions.  We  thus  apply  them 
in  a  loop  until  no  offending  triangles  remain: 

void  VertexInsertion: 
while  3  triangle  abc  violating  [U]  do 
ADD(x,  abc) 
endwhile. 

The  details  of  the  algorithm  for  adding  x  are  discussed  below. 

Consider  next  an  edge  ab  that  gets  too  short,  that  is,  Rab  =  ( C/Q)gab  at  time  tj.  The 
moment  in  time  tj  depends  on  points  a ,  b  and  their  length  scales  g(a ),  g(b),  which  all 
change  continuously  with  time.  To  remedy  the  violation  of  the  Lower  Size  Bound,  we 
contract  ab  by  removing  the  vertex  b  with  larger  length  scale  from  the  triangulation.  The 
removal  of  b  may  possibly  create  new  edges  violating  [L],  and  it  can  certainly  create 
triangles  violating  [U].  We  repair  the  triangulation  in  two  nested  loops. 

void  EDGECONTRACTION: 

while  3  edge  ab  violating  [L]  do 
if  g(a)  >  g(b)  then  a  ++  b  endif ; 

Remove  (£);  VertexInsertion 
endwhile. 


We  note  that  edge  contractions  and  vertex  insertions  can  also  be  used  to  modify  the 
restricted  Delaunay  triangulation  of  an  oversampling  until  it  satisfies  both  Size  Bounds. 
This  is  because  eveiy  oversampling  is  also  an  ^-sampling  and  thus  the  General  Home- 
omorphism  Theorem  applies.  However,  if  we  start  with  an  undersampling,  then 
the  algorithm  may  fail  because  conditions  needed  for  its  correct  operation  can  be 
violated. 

Vertex  Insertion.  Let  L  be  the  line  of  points  at  equal  distance  from  a,  b ,  c.  It  intersects 
the  plane  of  abc  in  the  circumcenter  z  of  the  triangle,  and  it  intersects  F  in  two  or  more 
points  of  which  we  add  the  point  x  e  L  D  F  closest  to  z.  We  prove  in  Section  8  that 
the  distance  between  x  and  z  is  necessarily  small  compared  with  gabc.  However,  even 
though  x  and  z  are  rather  close,  they  may  lie  in  different  mixed  cells.  To  find  jc,  we  start 
at  z,  determining  the  mixed  cell  that  contains  it  by  walking  from  a.  Then  we  walk  along 
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L  to  at  most  the  distance  specified  in  the  Circumcenter  Lemma  in  both  directions  from 
z.  Each  step  in  the  walk  enters  a  new  mixed  cell  fi.  Let  y  be  the  point  closest  to  z  that 
belongs  both  to  L  and  the  sphere  or  hyperboloid  of  p,.  If  y  lies  outside  (i  the  search 
continues.  Otherwise,  jc  =  y  and  we  add  x  to  the  triangulation  by  connecting  it  to  the 
edges  of  abc. 

Adding  x  to  D  as  described  is  likely  to  compromise  the  dual  correspondence  to 
the  restricted  Voronoi  diagram.  Inspired  by  the  incremental  algorithm  for  constructing 
Delaunay  triangulations  in  R2  [13],  we  repair  the  correspondence  by  edge  flipping.  More 
specifically,  we  push  each  edge  in  the  link  of  jc  on  an  empty  stack  and  then  process  the 
stack  until  it  runs  empty.  Let  pq  be  the  top  edge  on  the  stack  shared  by  triangles  xpq 
inside  and  pqy  outside  the  star  of  x.  Depending  on  a  local  geometry  test,  we  either  leave 
pq  as  an  edge  in  the  triangulation  or  we  flip  it  as  described  in  Section  5.  In  the  latter 
case  we  push  the  new  link  edges  py ,  yq  on  the  stack. 

void  ADD(x,  abc): 

assert  INSPHERE(*,  abc)\ 

substitute  xab ,  xb,  xbcy  xc,  xca ,  xa,  x  for  abc ; 

Push3  {ab,  be ,  ca ); 
while  not  ISEMPTY  do 
pq  =  POP; 

if  INSphereO,  pqy)  then 
EdgeFlipO#);  PusH2(/?y,  yq) 

endif 

endwhile. 

Let  z  be  the  dual  vertex  of  pqy  in  D  before  x  was  inserted.  Function  inSphere  decides 
whether  or  not  x  lies  inside  the  sphere  with  center  z  that  passes  through  p,q,y-  If  it 
does,  then  pqy  loses  the  reason  for  its  existence  which  justifies  the  flip. 

Edge  Contraction.  We  contract  an  edge  ab  by  removing  b  from  the  triangulation,  as 
illustrated  in  Fig.  10.  The  operation  removes  b  together  with  all  edges  and  triangles  in 
its  star,  and  it  covers  the  thus  created  hole  by  a  triangulation  without  interior  vertices. 
The  boundary  of  the  hole  is  a  topological  circle,  which  we  refer  to  as  a  polygon.  For 
each  vertex  p  let  p~  and  p+  be  its  predecessor  and  successor  in  an  ordering  around  the 
polygon.  The  algorithm  converts  the  star  of  b  into  the  new  triangulation  by  creating  a 
triangle  p~pp+  at  a  time  by  flipping.  When  only  three  triangles  remain  in  the  star  of  b 


Fig.  10.  The  removal  of  b  replaces  the  star  of  b  by  the  dotted  polygon  triangulation. 
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we  replace  its  star  by  a  single  triangle.  To  get  started,  we  push  all  vertices  in  the  link  of 
b  on  an  empty  stack. 

void  REMOVE(fc): 

for  all  p  in  link  of  b  do  PUSH (p)  endf or; 
while  stack  contains  more  than  three  vertices  do 
p  =  Pop; 

if  IND(/?“/?/?+)  then 
EdgeFlip(£>/?);  drop  p  from  polygon; 
if  p~  not  on  stack  then  PuSH(p“)  endif ; 
if  not  on  stack  then  PUSH (p+)  endif 

endif 
endwhile; 
p,q,r  =  Pop3; 

substitute  pqr  for  bpq,  bq ,  bqr ,  hr ,  brp ,  bp ,  b. 

Function  inD  returns  true  iff  all  other  vertices  of  the  current  polygon  lie  outside  the 
circumsphere  of  p ,  p+  whose  center  is  the  dual  restricted  Voronoi  vertex. 


7.  Topology  Adaptation 

The  way  the  skin  surface  is  connected  can  change  during  deformation.  This  section 
studies  when,  where,  and  how  these  changes  happen  in  the  growth  model.  It  also  de¬ 
scribes  how  we  locally  modify  the  general  sampling  strategy  to  avoid  the  computational 
impossibility  of  sampling  infinitely  many  points  accumulating  at  locations  of  infinite 
curvature. 

Growth  Model.  We  recall  that  the  growth  model  of  deformation  is  defined  by  changing 
the  square  radius  of  a  sphere  (a,  A)  from  A2  at  time  0  to  A2  -ft  at  time  t  e  R.  Compu¬ 
tationally,  this  is  the  simplest  kind  of  deformation  because  it  keeps  the  mixed  complex 
invariant.  Each  mixed  cell  contains  a  possibly  empty  sphere  or  hyperboloid  patch  of  the 
skin  surface.  After  normalization,  the  equation  of  the  sphere  or  hyperboloid  at  time  t  is 

x2  +  x\  ±  x\  =  ±R2  + 

Compare  this  with  (1)  and  (2)  in  Section  2.  A  metamorphosis  happens  when  the  right- 
hand  side  vanishes  at  time  t  =  q p2R2,  and  it  happens  at  the  center  but  only  if  the  center 
lies  in  the  interior  of  its  mixed  cell.  If  the  center  lies  outside,  the  portion  of  the  sphere 
or  hyperboloid  that  passes  through  the  center  is  not  part  of  and  thus  does  not  affect  the 
skin  surface.  The  special  case  where  the  center  lies  on  the  boundary  of  its  mixed  cell  is 
interesting.  We  will  see  that  in  this  case  the  metamorphosis  does  not  happen,  but  we  still 
have  to  modify  the  sampling  strategy  because  the  curvature  grows  beyond  any  bound. 

Using  local  considerations,  we  can  reduce  the  list  of  metamorphoses  to  the  four  given 
in  Table  1,  Section  2.  Cases  k  =  0,  3  correspond  to  an  appearing/disappearing  sphere. 
Cases  k  =  1,2  correspond  to  switching  a  hyperboloid  from  two  sheets  to  one,  or  vice 
versa.  In  each  case  we  can  interpret  the  center  as  a  critical  point  of  the  map  M:  R3  — ►  R 
whose  level  sets  M~l(t)  are  the  skin  surfaces  at  time  t.  Cases  k  =  0,  3  correspond  to 
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minima  and  maxima,  and  Cases  k  =  1,  2  to  two  types  of  saddle  points.  The  gradient 
of  M  vanishes  at  all  these  points  and  also  at  centers  that  lie  on  the  boundary  of  their 
mixed  cells.  The  latter  centers  correspond  to  degenerate  critical  points  in  the  sense  that 
an  arbitrarily  small  perturbation  of  M  suffices  to  turn  them  into  regular  points. 

Hot  Spots.  Common  to  every  metamorphosis  is  the  local  drop  in  length  scale,  which 
reaches  zero  at  the  moment  and  point  of  the  metamorphosis.  We  analyze  the  situation  in 
some  detail.  Let  H  be  a  positive  real  number.  The  hot  portion  of  the  skin  surface  F  is 
the  set  of  points  with  length  scale  H  or  smaller, 

Fh  =  {*  €  F  |  e(x)  <  H}. 

By  the  Iso-curvature  Lemma,  we  have  g(x)  <  H  only  if  x  is  sufficiently  close  to  the 
center  of  a  sphere  or  hyperboloid.  Let  zx  be  such  a  center.  We  call  the  ball  fix  =  {y  £ 
R3  I  \\y  —  Zx\\  <  H}  the  hot  ball  of  X .  A  hot  ball  is  relevant  only  inside  its  mixed  cell. 
The  union  of  hot  balls,  each  clipped  to  within  its  mixed  cell,  is  the  hot  portion  of  space, 
denoted  as  R^. 

Hot  Spot  Lemma.  Fh  =  F  C 1  R^. 

In  words,  a  point  x  e  F  belongs  to  the  hot  portion  of  the  skin  surface  iff  it  belongs  to 
the  hot  portion  of  space.  In  the  growth  model  the  hot  portion  of  space  is  constant,  while 
the  hot  portion  of  the  skin  changes  as  the  surface  moves  through  that  portion  of  space. 
The  Hot  Spot  Lemma  follows  directly  from  the  Iso-curvature  Lemma  and  does  not  need 
a  separate  proof. 

Depending  on  whether  centers  lie  inside  or  outside  their  mixed  cells,  the  hot  portion  of 
space  is  locally  a  union  or  intersection  of  hot  balls.  The  mixed  complex  decomposes  this 
union  and  intersection  into  convex  pieces,  as  illustrated  in  Fig.  11.  The  common  radius 
of  all  hot  balls  is  H.  As  long  as  none  of  the  centers  lies  on  the  boundary  of  its  mixed 
cell,  we  can  eliminate  any  overlap  by  decreasing  H  while  keeping  it  positive.  We  will 


Fig.  11.  Dotted  Voronoi  diagram,  dashed  Delaunay  triangulation,  solid  mixed  complex,  solid  data  points, 
hollow  other  centers,  and  shaded  hot  portion  of  space.  Each  label  shows  the  dimension  of  the  Delaunay  simplex 
involved  in  the  construction  of  the  mixed  cell. 
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shortly  discover  that  an  even  stronger  separation  property  between  hot  balls  is  needed 
to  prevent  edges  reaching  from  one  to  another,  which  can  be  achieved,  e.g.,  by  choosing 
H  equal  to  half  the  value  that  guarantees  pairwise  disjointness.  A  center  on  a  mixed  cell 
boundary  has  probability  zero  and  is  considered  a  degenerate  case.  For  now,  we  simplify 
the  discussion  by  assuming  the  non-degenerate  case,  where  H  is  small  enough  such  that 
hot  balls  are  pairwise  disjoint.  We  will  return  to  the  degenerate  case  shortly. 

Time  for  Change.  The  hot  portion  is  more  difficult  to  triangulate  than  the  rest  of  the 
skin  surface.  One  reason  is  the  metamorphosis,  another  is  the  accumulation  of  vertices 
in  a  small  region.  The  sphere  case  is  relatively  harmless,  because  the  area  decreases  at 
the  same  rate  as  the  density  requirement  increases.  Indeed,  a  constant  number  of  vertices 
suffices  to  shrink  a  sphere  to  an  arbitrarily  small  size.  The  case  of  a  hyperboloid  that 
approaches  its  limiting  double-cone  is  more  problematic,  because  the  number  of  vertices 
near  the  center  grows  beyond  any  bound.  To  circumvent  the  computational  impossibility 
of  sampling  infinitely  many  points,  we  change  the  sampling  strategy  inside  the  hot 
balls.  We  give  up  on  e-sampling  to  get  a  sparse  sampling,  but  we  preserve  the  closed 
ball  property.  The  triangulation  algorithm  remains  oblivious  to  the  changed  sampling 
density  and  keeps  constructing  the  restricted  Delaunay  triangulation. 

Consider  a  two-sheeted  hyperboloid  and  translate  time  such  that  the  metamorpho¬ 
sis  happens  at  time  t  =  0.  The  hyperboloid  enters  its  hot  ball  at  time  — 2//2,  turns 
into  a  double-cone  at  time  0,  and  leaves  the  hot  ball  as  a  one-sheeted  hyperboloid  at 
time  2 H2.  The  special  sampling  strategy  that  allows  us  to  go  through  this  motion  de¬ 
pends  on  a  parameter  0  <  h  <  1.  Special  sampling  begins  at  time  /0  =  —2H2h2 
when  the  two-sheeted  hyperboloid  enters  the  ball  of  radius  Hh,  and  it  ends  at  time 
t\  —  2H2h2  when  the  one-sheeted  hyperboloid  leaves  that  ball,  as  shown  in  Fig.  12. 
At  time  fy,  the  hyperboloid  intersects  the  boundary  of  the  hot  ball  in  two  hot  cir¬ 
cles.  The  shape  adaptation  algorithm  moves  these  circles  along  their  integral  lines, 
which  implies  that  they  grow  from  radius  R0  =  Hy/(\  -  h2)/2  at  time  f0  to  radius 
Ri  =  Hyf(\  +  h2)/2  at  time  t\.  Simultaneously,  the  distance  between  the  two  cir¬ 
cles  decreases  from  2R\  to  2R o.  We  define  the  hot  sphere  to  pass  through  the  two 
hot  circles.  At  time  to  and  t\ ,  it  is  the  boundary  of  the  hot  ball,  but  in  the  open  time 
interval  between  t0  and  t\ ,  it  is  cocentric  and  smaller  than  that  boundary.  General  sam¬ 
pling  applies  outside  the  hot  sphere  and  special  sampling  applies  on  and  inside  that 
sphere. 


Fig.  12.  Head-on  view  of  start,  middle,  end  configurations  generated  by  special  sampling  taking  a  two- 
sheeted  to  a  one-sheeted  hyperboloid.  The  hot  sphere  is  solid  and  the  sphere  that  triggers  the  metamorphosis 
is  dotted. 
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Extreme  Configuration .  At  time  to,  we  kick  off  special  sampling  by  creating  the  double¬ 
cup  as  the  start  configuration  of  the  metamorphosis  representing  the  intruding  portion 
of  the  two-sheeted  hyperboloid,  as  shown  in  the  left  drawing  of  Fig.  12.  Consider  one 
sheet  of  the  hyperboloid  and  let  a  be  its  intersection  point  with  the  symmetry  axis.  Let 
bo,  b\, . . . ,  bt- 1  be  the  vertices  of  a  regular  £-gon  along  the  hot  circle  in  this  sheet. 
We  mirror  these  points  across  the  symmetry  plane  of  the  hyperboloid  and  get  points 
a! ,  b'Q,  b\ , . . .  on  the  other  sheet. 

Add2 (a,  af)\ 

for  /  =  0  to  £  —  1  do  ADD2(fy,  fcj)  endfor; 

EdgeContraction. 

Recall  that  Function  Add  really  takes  two  parameters,  namely  a  point  and  a  triangle 
whose  circumsphere  centered  at  the  dual  restricted  Voronoi  vertex  encloses  the  point. 
Each  call  to  the  function  must  therefore  be  preceded  by  a  search  for  such  a  triangle. 
Whenever  we  contract  an  edge  by  removing  one  of  its  endpoints  we  make  sure  that  the 
endpoint  is  not  one  of  the  newly  added  vertices.  Section  10  derives  sufficient  conditions 
for  h  and  £  that  guarantee  the  above  algorithm  successfully  constructs  the  double-cup  as 
the  start  configuration  of  the  metamorphosis.  By  this  we  mean  that 

(i)  a,  a'  are  the  only  vertices  inside  and  bt,  b\,  for  0  <  i  <  £,  are  the  only  vertices 
on  the  hot  sphere, 

(ii)  the  link  of  a  in  D  is  the  regular  f-gon  of  vertices  bi9  and  symmetrically  the  link 
of  af  is  the  f-gon  of  vertices  b\ . 

Assuming  C  =  0.08,  Q  =  1.65  we  will  see  that  h  =  0.98  and  £  =  5  are  feasible  values 
for  the  two  constants.  For  ease  of  reference  we  say  the  vertices  and  edges  in  the  links  of 
a ,  a'  are  hot  and  the  vertices,  edges,  and  triangles  in  the  stars  of  a,  a'  are  very  hot . 

The  end  configuration  of  the  metamorphosis  is  similar  to  the  start  configuration  of 
the  inverse  metamorphosis.  As  shown  in  the  right  drawing  of  Fig.  12,  it  consists  of 
two  rings  of  triangles  forming  a  cylinder-with-a-waist  representing  the  intruding  portion 
of  the  one-sheeted  hyperboloid.  Let  u0,  u\,  •  •  -  ■  um-\  be  the  vertices  of  a  regular  m- 
gon  along  the  waist  where  the  hyperboloid  intersects  its  symmetry  plane.  Similarly,  let 
Wo,  W\ , . . . ,  wm-i  be  the  vertices  of  another  regular  m-gon  along  one  of  the  two  hot  cir¬ 
cles,  rotated  by  n /m  relative  to  the  m-gon  along  the  waist.  Finally,  let  wf0,  w[ , . . . , 
be  the  vertices  of  the  mirror  m-gon  on  the  other  hot  circle. 

for  j  =  0  to  m  —  1  do  ADD3(Mj,  Wj,  wf.)  endfor; 

EdgeContraction. 

As  before  we  search  for  an  offending  triangle  before  we  add  a  point,  and  we  make  sure 
that  the  removed  vertex  of  every  contracted  edge  pq  is  not  one  of  the  newly  added 
ones.  Section  10  derives  sufficient  conditions  for  h  and  m  that  guarantee  the  above 
algorithm  successfully  constructs  the  cylinder-with-a-waist  as  the  start  configuration  of 
the  inverse  metamorphosis.  What  precisely  we  mean  by  this  should  be  obvious.  Assuming 
C  =  0.08,  Q  =  1.65  we  will  see  that  h  =  0.98  and  m  —  40  are  feasible  values  for  the 
two  constants.  For  ease  of  reference  we  again  say  that  the  vertices  and  edges  along  the 
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two  hot  circles  are  hot  and  that  the  vertices,  edges,  and  triangles  between  the  two  hot 
circles  are  very  hot. 

In  the  forward  direction  we  switch  from  the  double-cup  to  the  cylinder-with-a-waist 
at  time  0,  and  in  the  backward  direction  we  do  it  the  other  way  round.  The  latter  is  easier 
because  we  just  need  to  meld  the  m-gon  of  the  waist  into  a  single  vertex  and  then  split 
that  vertex  into  two.  In  the  forward  direction  we  first  meld  a  and  af  into  a  single  vertex 
and  then  expand  that  vertex  into  a  regular  polygon  (interleaving  angularly  between  the 
b{).  The  expansion  creates  two  new  rings  of  triangles  between  the  new  polygon  and  the 
polygons  representing  the  two  hot  circles.  This  is  done  following  the  angular  order  of 
the  involved  vertices  around  the  symmetry  axis. 

Special  Sampling.  The  main  difference  between  special  and  general  sampling  is  that 
the  former  gives  up  on  the  Lower  Size  Bound  for  hot  edges  and  on  the  Upper  Size  Bound 
for  very  hot  triangles.  The  length  of  hot  edges  is  bounded  from  above  because  the  Upper 
Size  Bound  applies  to  the  incident  triangles  outside  the  hot  sphere.  A  more  detailed 
analysis  of  edge  and  triangle  sizes  including  a  proof  of  the  closed  ball  property  in  spite 
of  special  sampling  is  given  in  Section  10. 

The  goal  of  special  sampling  is  to  maintain  the  double-cup  and  the  cylinder-with- 
a-waist  during  the  first  and  the  second  halves  of  the  time  interval.  It  acts  primarily 
by  modifying  general  sampling  for  points  on  and  inside  the  hot  sphere.  As  a  general 
rule,  an  edge  is  contracted  by  removing  an  endpoint  that  is  not  hot.  Cases  where  both 
endpoints  are  hot  occur  only  at  the  end  of  the  metamorphosis  (or  its  inverse)  and  will 
be  discussed  separately.  There  are  two  ways  in  which  general  sampling  can  intrude  into 
the  hot  sphere:  by  adding  a  point  inside  that  sphere  and  by  flipping  a  hot  edge.  In  both 
cases  we  prevent  the  intrusion  by  bisecting  the  endangered  hot  edge  be ,  as  illustrated  in 
Fig.  13.  Specifically,  we  add  the  midpoint  q  of  the  shorter  hot  circle  arc  that  connects  b 
with  c.  The  addition  of  q  may  create  edges  that  violate  the  Lower  Size  Bound.  Of  these 
we  contract  the  ones  that  are  not  hot,  always  making  sure  we  remove  the  endpoint  that  is 
not  hot.  As  discussed  above,  we  choose  H  small  enough  so  that  hot  spheres  cannot  get 
too  close  to  each  other  and  every  non-hot  edge  has  at  least  one  non-hot  vertex.  Infinite 
loops  cannot  occur  because  each  iteration  leaves  an  additional  hot  vertex  behind.  The 
hot  circle  gets  denser  and  intrusions  into  the  hot  sphere  get  progressively  more  difficult. 

Special  sampling  maintains  the  special  configuration,  but  it  does  not  guarantee  the 
two  Size  Bounds.  They  must  therefore  be  enforced  algorithmically  at  the  end  of  the 
metamorphosis. 


SpecialVertexInsertion; 
Edge  Contraction. 


Fig.  13.  The  hot  edge  be  is  bisected  either  because  the  dual  restricted  Voronoi  vertex  of  bed  lies  inside  the 
hot  sphere  or  edge  flipping  attempts  to  change  be  to  da. 
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The  clean-up  operation  is  correct  if  we  maintain  the  closed  ball  property,  which  is  initially 
guaranteed  by  special  sampling.  While  maintaining  that  property  might  be  difficult  in 
general,  we  can  use  the  insights  gained  from  the  proofs  of  the  two  Persistence  Lemmas 
in  Section  10  and  add  points  only  on  the  two  hot  circles.  This  is  the  difference  between 
the  above  function  and  Function  VERTEXlNSERTlON  introduced  earlier.  The  size  analysis 
in  Section  10  implies  that  we  can  satisfy  the  Upper  Size  Bound  even  with  this  restriction 
on  new  vertex  locations. 

Degenerate  Centers.  Recall  that  a  degenerate  center  is  one  that  lies  on  the  boundary 
of  its  mixed  cell.  Each  facet  lies  half-way  between  the  centers  of  the  two  mixed  cells 
that  share  it.  If  it  contains  one  center,  then  it  also  contains  the  other,  which  implies 
that  a  degenerate  center  is  also  a  multiple  center.  In  a  non-degenerate  mixed  complex 
(which  may  still  have  degenerate  centers)  every  facet  is  shared  by  two,  every  edge  by 
four,  and  every  vertex  by  eight  mixed  cells.  In  the  Morse  theoretic  view  of  centers 
as  critical  points,  each  degenerate  center  is  the  location  where  several  critical  points 
collide  and  cancel  each  other.  This  is  why  there  is  no  metamorphosis  at  those  loca¬ 
tions.  When  the  surface  moves  through  such  a  degenerate  center,  its  curvature  does 
blow  up  momentarily,  but  the  surface  then  becomes  smooth  again,  with  no  topology 
change. 

The  remainder  of  this  section  describes  the  various  types  of  degenerate  centers.  The 
enumeration  exhausts  the  cases  that  occur  in  non-degenerate  mixed  complexes.  Consis¬ 
tent  with  Table  1,  the  label  of  a  mixed  cell  is  k  if  it  is  constructed  from  a  ^-dimensional 
Delaunay  simplex  and  its  dual  (3  —  k) -dimensional  Voronoi  polyhedron.  We  label  the 
facets,  edges,  and  vertices  by  concatenating  the  labels  of  the  mixed  cells  that  share  them. 
There  are  three  facet  types  labeled  01,  12,  23,  two  edge  types  labeled  01 12,  1223,  and 
one  vertex  type  labeled  01 112223. 

The  case  where  the  degenerate  center  lies  in  the  interior  of  a  01  facet  is  illustrated 
in  Fig.  14.  Reading  the  top  row  of  the  figure  from  left  to  right  we  see  the  skin  surface 
passing  through  the  center.  The  bottom  row  shows  the  hot  ball  around  the  degenerate 
center  and  its  intersection  with  the  evolving  surface  and  the  body  it  bounds.  The  case  of  a 
degenerate  center  on  a  23  facet  is  symmetric.  The  shrunken  Voronoi  polygon  is  replaced 


Fig.  14.  Degenerate  center  in  the  interior  of  a  01  facet.  Evolution  of  skin  surface  at  the  top  and  of  hot  ball  at 
the  bottom.  The  shaded  regions  on  the  hot  spheres  show  the  intersection  with  the  body. 
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Fig.  15.  Degenerate  center  in  the  interior  of  a  1 2  facet.  Evolution  of  skin  surface  at  the  top  and  of  hot  sphere 
at  the  bottom.  The  hyperboloid  in  the  type-2  cell  has  a  vertical  symmetry  axis,  while  the  hyperboloid  in  the 
type-1  cell  behind  it  has  a  horizontal  one. 


by  a  shrunken  Delaunay  triangle,  and  we  read  Fig.  14  from  right  to  left.  Furthermore, 
the  body  lies  on  the  other  side  of  the  skin  surface. 

The  case  of  a  degenerate  center  in  the  interior  of  a  12  facet  is  different  and  geometri¬ 
cally  more  interesting.  Both  mixed  cells  contain  hyperboloids,  and  their  symmetry  axes 
pass  orthogonally  through  the  common  center.  Each  symmetry  axis  lies  in  the  symmetry 
plane  of  the  other  hyperboloid,  which  forms  a  right  angle  with  the  plane  of  the  facet.  The 
evolution  of  the  skin  surface  passing  through  such  a  degenerate  center  is  illustrated  in  the 
top  row  of  Fig.  15.  As  shown  in  the  bottom  row,  the  skin  surface  intersects  the  boundary 
of  the  hot  ball  around  the  degenerate  center  in  a  curve  consisting  of  four  half-circles. 
Instead  of  a  pair  of  hot  circles,  we  have  a  hot  curve  consisting  of  four  half-circles. 

The  case  of  a  degenerate  center  in  the  interior  of  an  edge  labeled  01 1 2  is  illustrated  in 
Fig.  16.  The  two  pairs  of  facets  common  to  a  type  1  cell  form  aright  dihedral  angle  each, 
and  the  remaining  two  dihedral  angles  add  up  to  n  but  are  otherwise  arbitrary.  As  shown 
to  the  right,  the  skin  surface  intersects  the  boundary  of  the  hot  ball  in  four  circular  arcs, 
two  of  which  are  half-circles.  The  angles  remain  the  same  during  the  entire  transition 


Fig.  16.  Degenerate  center  in  the  interior  of  a  01 12  edge.  Snapshot  of  skin  surface  and  hot  ball  at  the  moment 
the  surface  passes  through  the  center. 
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in  which  the  skin  surface  sweeps  through  the  hot  ball.  The  case  of  a  degenerate  center 
on  an  edge  labeled  1223  is  symmetric,  and  Fig.  16  applies,  except  that  the  body  lies  on 
the  other  side  of  the  surface.  In  the  case  of  a  degenerate  center  at  a  vertex,  the  hot  curve 
consists  of  three  pairs  of  circular  arcs.  The  angles  of  the  arcs  remain  the  same  during  the 
transition  in  which  the  skin  surface  sweeps  through  the  hot  ball. 


Part  HI.  Analysis 

The  next  three  sections  analyze  the  algorithm  and  the  triangulations  it  creates.  Section  8 
studies  questions  related  to  sampling  density,  Section  9  focuses  on  scheduling,  and 
Section  10  examines  the  topology  adaptation  algorithm. 


8.  Sampling  Density 

We  derive  conditions  for  the  constants  C  and  Q  in  order  to  prove  the  curvature  adaptation 
algorithm  in  Part  II  is  correct. 

Conditions.  We  prove  that  point  insertions  do  not  generate  edges  that  violate  the  Lower 
Size  Bound.  That  proof  requires  that  Q  is  not  too  large.  We  also  prove  that  the  restricted 
Voronoi  vertex  dual  to  a  triangle  can  be  found  near  the  circumcenter  of  that  triangle.  That 
proof  requires  that  the  vertices  of  the  triangulation  form  an  5-sampling.  Finally,  we  prove 
that  the  vertices  indeed  form  an  5-sampling,  with  5  satisfying  Condition  (I).  The  closed 
ball  property  established  in  Section  4  then  implies  that  the  triangulation  produced  is 
homeomorphic  to  the  skin  surface.  That  proof  relies  on  the  quality  of  the  approximation, 
which  is  guaranteed  by  the  algorithm  provided  C  Q  is  not  too  large.  For  ease  of  reference 
we  collect  the  conditions  before  deriving  them. 

(II)  Q2 -4CQ-2>0. 

(ID)  <$2/(l  4-  <$)2  -  84/ 4  >  C2Q2, 

where  <$  =  5-2C(5  +  l)/(2-i-  2C).  We  get  (II)  and  (III)  as  sufficient  conditions  for  the 
proofs  of  the  No-Short-Edge  Lemma  and  the  Sampling  Lemma  below.  Condition  (II) 
is  equivalent  to  Q  >  2C  +  \/4C2  +  2.  Assuming  5  =  50  =  0.279 . . . ,  we  can  satisfy 

Conditions  (II)  and  (HI)  by  setting  C  =  0.08  and  Q  =  1.65.  In  this  case  8  =  0.166 - 

Small  improvements  are  possible. 

Short  Edges.  An  edge  contraction  may  perhaps  cause  other  edge  contractions,  but  this 
cannot  go  on  forever  because  we  eventually  violate  the  Upper  Size  Bound.  Similarly,  a 
vertex  insertion  may  cause  other  vertex  insertions,  but  this  cannot  go  on  forever  because 
we  eventually  violate  the  Lower  Size  Bound.  It  is  possible  that  an  edge  contraction 
causes  vertex  intersections,  but  a  vertex  insertion  cannot  cause  edge  contractions.  This  is 
because  a  vertex  insertion  cannot  create  edges  of  size  below  the  allowed  threshold.  This 
is  what  prevents  infinite  loops  in  spite  of  the  algorithm’s  partially  conflicting  efforts  to 
avoid  short  edges  and  large  triangles  simultaneously.  Let  abc  be  the  triangle  that  causes 
the  addition  of  the  dual  restricted  Voronoi  vertex  x  e  F. 
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No-Short-Edge  Lemma.  Every  edge  xy  created  during  the  addition  of  x  has  size 
larger  than  (C/Q)gxy. 


Proof  We  have  Rabc  >  CQQabc ■  The  sphere  with  center  x  that  passes  through  a,  b , 
c  has  radius  X  >  Rabc  and  it  contains  no  vertex  other  than  x  inside.  Every  new  edge 
xy  has  therefore  length  \\x  —  y||  >  X  >  CQgabc •  Assume  without  loss  of  generality 
that  Qabc  =  Q (a).  We  use  the  Curvature  Variation  Lemma  to  derive  upper  bounds  for  the 
length  scales  at  x  and  y: 


QW  <  Q(a)  +  X  <  ^-2-  +  ||*  _  y\\t 

e(y)  <  e(x)  +  \\x  -  y||  <  ^-2-  +  2^  \\x  -  y\\. 


Hence 


RXy  — 


\\x-y\\  max(e(jc),  e(y)} 


Condition  (II)  implies  C/Q 
claimed. 


2  "  4  +  2/C0  * 

CQ/(4CQ  -f  2),  and  therefore  Rxy 


( C/Q)Qxy ,  as 
□ 


Close  Dual  Vertices.  Consider  the  point  addition  triggered  by  the  triangle  abc  violating 
the  Upper  Size  Bound.  As  before,  we  denote  the  line  of  points  at  equal  distance  from  a, 
b ,  c  by  L,  the  circumcenter  of  abc  by  z,  and  the  point  ofLHF  closest  to  z  by  x.  We 
prove  an  upper  bound  on  the  distance  between  x  and  z  assuming  an  ^-sampling  of  F. 

Circumcenter  Lemma.  The  distance  between  x  and  z  is  \\x  —  z||  <  (e2/2 )Qabc- 


Proof  Assume  Qabc  =  g{a)  <  e(b),Q(c).  We  have  \\x  -  a ||  <  £q(x)  by  assump¬ 
tion  of  ^-sampling  and  therefore  (1/(1  4-  e))Qabc  <  £?(*)  by  the  Curvature  Varia¬ 
tion  Lemma.  We  get  an  upper  bound  on  the  distance  between  x  and  z  by  assuming 
f?(*)  is  as  small  as  possible  and  a,  b,  c  lie  on  the  sandwiching  sphere  with  radius 
Q(x)  =  (1/(1  -1-  e))Qabc  passing  through  x .  This  configuration  is  sketched  in  Fig.  17. 
Note  that  ||*  —  z\\/\\x  —  a\\  =  ||*  —  a\\/2g(x)  by  equality  of  angles  formed  by  orthog¬ 
onal  sides.  Therefore, 


as  claimed. 


II*- 


ii*  -an2 
2 <?(*) 


<  jeM. 


□ 


Fig.  17.  Dashed  sphere  of  radius  q(x )  passing  through  a,b,cyx  and  bold  circle  with  center  z  passing  through 
a,  by  c. 
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The  relevance  of  the  Circumcenter  Lemma  to  the  curvature  adaptation  algorithm 
should  be  obvious.  When  the  triangle  abc  violates  Condition  (III),  we  need  first  to  find 
its  dual  vertex  in  the  restricted  Voronoi  diagram  and  then  add  this  vertex  to  V.  This 
vertex  is  the  point  x,  and  the  Circumcenter  Lemma  gives  a  bound  on  how  far  from  z  we 
have  to  search  before  we  are  guaranteed  to  find  x.  As  shown  in  the  proof  of  the  Voronoi 
Edge  Lemma,  each  additional  point  y  e  L  n  F  is  too  far  from  z  to  possibly  belong  to 
the  Voronoi  edge  dual  to  abc. 

Maintaining  Density.  We  show  that  the  algorithm  for  curvature  adaptation  maintains 
the  ^-sampling  property  of  the  vertex  set.  Recall  that  this  means  that  for  every  point 
x  e  F  there  is  a  vertex  a  e  V  whose  distance  from  x  is  ||a  —  x  ||  <  e £>(*).  The  constant 
e  is  to  be  chosen  so  it  satisfies  Condition  (I). 

It  is  interesting  to  see  that  the  two  Size  Bounds  by  themselves  are  too  weak  to  imply 
^-sampling.  We  can  put  four  points  near  each  other  on  a  sphere  in  such  a  way  that 
all  four  triangles  and  six  edges  satisfy  [L]  and  [U].  Nevertheless,  the  boundary  of  the 
tetrahedron  is  a  miserably  inadequate  approximation  of  the  sphere  surface.  We  argue 
that  the  algorithm  cannot  get  to  this  problematic  state,  because  of  the  way  it  would 
temporarily  have  to  violate  the  two  Size  Bounds.  In  other  words,  we  use  continuity  in 
time  to  prove  the  claim  on  sampling.  In  stating  the  result,  we  assume  the  skin  surface 
deforms  continuously  with  time.  For  now  we  disallow  metamorphoses.  Let  t0  <  t\  be 
two  points  in  time  so  the  topological  type  is  constant  within  |>o,  h].  We  write  F(t )  for 
the  skin  surface  at  time  t  and  Vo,  Vi  for  the  vertex  sets  at  times  t0,  t\. 

Sampling  Lemma.  If  V0  is  an  E-sampling  ofF(to ),  then  V\  is  an  E-sampling  ofF(t\). 


Proof.  Assume  the  opposite  and  let  t  e  [to,h]  be  the  first  moment  in  time  when  the 
skin  surface  is  not  e-sampled.  Then  there  is  a  point  x  e  F(t)  such  that  no  vertex  lies 
inside  the  sphere  with  center  x  and  radius  X  =  eq(x).  By  minimality  of  t,  the  sphere 
passes  through  at  least  one  vertex,  a ,  but  we  need  three.  To  get  two  more,  we  continuously 
increase  the  sphere  while  keeping  its  center  on  the  surface.  Vertex  a  remains  on  the  sphere 
at  all  times  and  we  permit  no  vertices  inside  the  sphere.  Let  y  e  F(t)  be  the  center  when 
we  reach  the  other  two  vertices,  b  and  c.  The  radius  of  the  new  sphere  is  Y  >  X  because 
the  radius  can  only  increase  from  x  to  y.  Using  the  Curvature  Variation  Lemma,  we  get 
Qabc  <  £(x)(l  +  e)  and  therefore  (s/(  1  +  e))Qabc  <  eq(x)  <  Y.  Assume  without  loss 
of  generality  that  Qabc  —  Q(a),  and  let  z  be  the  circumcenter  of  abc.  The  Upper  Size 
Bound  implies  ||z  —  a\\  =  Rabc  <  C  Qg(a).  Using  the  Circumcenter  Lemma,  we  get  an 
upper  bound  on  the  square  distance  between  y  and  a , 

Y 2  =  lly  -zll2  +  Ik  -  all2 

<  jeL  +  c2e2eL- 


This  implies 


7rrv<T  +  c2Q2’ 
(1  +  e)2  4 


which  contradicts  Condition  (III). 


□ 
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For  example  for  C  =  0.08,  Q  =  1.65  the  Sampling  Lemma  holds  for  all  e  in  an 

interval  with  endpoints  0.15 .. .  and  0.98 _ 

The  algorithm  explained  in  Section  7  maintains  an  ^-sampling  across  metamorphoses. 
More  precisely,  it  violates  the  required  sampling  density  inside  the  hot  sphere  of  each 
metamorphosis.  As  proved  in  Section  10,  the  special  sampling  strategy  repairs  the  e- 
sampling  property  before  the  skin  surface  comes  out  of  the  hot  sphere,  and  it  maintains 
the  closed  ball  property  at  all  times.  The  Sampling  Lemma  thus  generalizes  to  any  time 
interval,  including  those  that  contain  metamorphoses. 

Denser  Sampling.  Even  though  the  Sampling  Lemma  proves  that  the  vertices  form  an 
e-sampling  between  all  operations,  it  allows  for  momentary  violations  of  the  density 
requirement  during  the  contraction  of  an  edge  ab.  Specifically,  e-samplings  are  not 
guaranteed  right  after  b  is  removed  and  before  appropriate  vertex  insertions  repair  the 
Upper  Size  Bound.  The  contraction  happens  only  if  the  Lower  Size  Bound  is  violated, 
which  implies  || b  -  a\\  <  (2 C/(Q)Qab.  Since  such  an  edge  contraction  always  removes 
the  vertex  with  larger  length  scale  we  have  g(b)  >  g(a)  and  therefore  Qab  <  e(x)  + 
Ik  —  b\\.  For  every  point  x  €  F  whose  only  vertex  within  distance  eq(x)  is  b  we  thus  have 

11^ -all  <  Ik  -b\\  +  \\b-a\\ 

2  C 

<  €Q(X)  +  —  (Q(X)  +  €Q(X)) 

2  C  1 

=  £  +  ~q^£  +  i)J  e(*)« 

In  words,  the  vertices  still  form  an  e -sampling,  but  for  a  somewhat  larger  value  of  e.  To 
say  the  same  thing  in  reverse  we  define  S  =  e  -  ( 2C/(Q  +  2 C))(e  +  1)  and  note  that 
£  =  S  4-  (2C/Q)(8  +  1).  If  we  choose  the  constants  C  and  Q  such  that  the  Sampling 
Lemma  implies  the  maintenance  of  a  <5 -sampling,  then  we  have  an  ^-sampling  at  all 
times,  even  during  the  execution  of  an  edge  contraction.  Condition  (III)  enforces  such  a 
choice  of  constants. 


9.  Scheduling 

The  overall  algorithm  deforms  the  skin  surface  by  executing  operations  ordered  in  time. 
Some  of  these  operations  require  others  to  repair  the  damage,  and  these  others  are 
executed  following  a  partial  rather  than  a  total  order.  As  a  general  rule,  total  ordering 
is  more  expensive  but  easier  to  prove  correct  than  partial  ordering.  This  section  reviews 
all  operations  and  discusses  their  treatment  by  the  scheduling  algorithm.  It  also  provides 
correctness  proofs  for  the  flipping  algorithms  used  to  restore  the  restricted  Delaunay 
triangulation  after  vertex  insertions  and  edge  contractions. 

Total  and  Partial  Ordering.  Operations  triggered  by  the  motion  of  the  skin  surface  are 
ordered  in  time.  We  have  five  types: 

1.  coordinate  updates, 

2.  edge  flips, 

3.  edge  contractions, 
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4.  vertex  insertions, 

5.  metamorphoses. 

Vertex  coordinates  change  continuously  with  time,  and  we  avoid  most  of  the  related 
computational  expense  by  updating  coordinates  when  and  only  when  they  are  used  by 
other  operations.  The  last  four  operations  are  discrete  events  that  are  stored  in  a  priority 
queue  ordered  by  time.  The  moment  in  time  when  an  edge  flip,  edge  contraction,  or  vertex 
insertion  matures  is  a  root  of  a  continuous  function.  In  the  growth  model  of  deformation, 
the  moment  in  time  when  a  metamorphosis  matures  is  predictable  from  the  ordering  of 
Delaunay  simplices  described  in  Section  2.  For  more  general  deformations,  the  time  of 
a  metamorphosis  is  also  a  root  of  a  continuous  function. 

Each  operation  other  than  the  coordinate  update  and  the  edge  flip  is  further  decom¬ 
posed  into  a  sequence  of  operations.  For  example  a  vertex  insertion  relies  on  point 
additions  and  edge  flips  to  achieve  the  desired  effect  locally  and  restore  the  restricted 
Delaunay  triangulation.  Conceptually,  such  a  sequence  is  executed  at  an  instant,  while 
time  stands  still.  We  cannot  therefore  resort  to  time  for  a  global  ordering  mechanism. 
The  operations  in  each  sequence  are  therefore  scheduled  following  a  partial  rather  than 
a  total  order.  The  most  frequently  executed  operation  is  the  edge  flip.  The  choice  of 
constants  C  and  Q  guarantees  that  the  restricted  Voronoi  diagram  has  the  closed  ball 
property  at  all  times,  even  in  the  middle  of  an  edge  contraction.  We  would  therefore 
expect  that  a  simple  iteration  of  edge  flips  will  suffice  to  restore  the  restricted  Delaunay 
triangulation.  While  this  is  easy  to  prove  for  point  additions,  it  is  possibly  incorrect  for 
point  removals.  This  is  why  we  resort  to  the  more  complicated  edge  flipping  algorithm 
for  point  removal  described  in  Section  6. 

Flipping  after  Point  Addition.  A  vertex  insertion  operation  is  recursive  and  unwinds 
into  a  sequence  of  point  additions  (usually  just  one),  each  followed  by  a  sequence  of 
edge  flips.  Let  D0  and  D\  be  the  restricted  Delaunay  triangulations  immediately  before 
and  after  adding  the  vertex  x.  We  focus  on  the  sequence  of  flips  following  the  addition 
of  x  and  argue  that  this  sequence  successfully  constructs  D\ . 

As  explained  in  Section  6,  the  algorithm  maintains  a  subset  of  the  edges  in  the  link 
of  x  on  a  stack.  Each  edge  pq  in  the  link  belongs  to  two  triangles,  xpq  inside  and  pqy 
outside  the  star  of  x.  Thus  pqy  exists  in  D0  and  it  remains  in  D\  iff  x  lies  outside  the 
circumscribed  sphere  of  pqy  whose  center  is  the  dual  restricted  Voronoi  vertex.  The 
INSPHERE  test  used  to  decide  whether  or  not  to  flip  pq  captures  exactly  that  information. 
If  the  decision  is  negative  we  know  that  pqy  and  pq  remain  in  D\ .  Otherwise,  the  flip  of 
pq  increases  the  star  of  x  by  one  triangle,  and  it  decreases  the  portion  of  the  triangulation 
outside  the  star  by  one  triangle.  The  monotonicity  of  the  transfer  implies  that  the  loop 
of  edge  flips  halts.  There  are  no  obstacles  to  flipping  pq  other  than  if  xy  is  already  an 
edge  of  the  triangulation,  which  is  the  case  iff  either  p  or  q  belongs  to  only  three  edges. 
Say  p  belongs  to  px,  pq,  py  and  to  no  other  edges.  By  the  closed  ball  property  of  the 
restricted  Voronoi  diagram,  p  belongs  to  at  least  three  edges  of  D\,  and  because  p  can 
only  lose  edges  during  flipping,  the  inSphere  test  for  x,  pqy  must  be  negative  and  the 
flip  of  pq  will  not  be  attempted.  Another  conceivable  reason  for  failure  is  that  the  link  of 
x  does  not  reach  a  triangle  that  should  be  removed.  However,  this  is  impossible  because 
the  closed  ball  property  implies  that  Dq  —  D\  is  an  open  disk,  and  the  link  of  x,  which  is  a 
topological  circle,  is  adj  acent  to  triangles  in  Z)0  —  A  until  it  has  swept  out  that  entire  disk. 
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Flipping  after  Point  Removal  An  edge  contraction  unwinds  into  a  sequence  of  point 
removals  and  point  additions,  each  followed  by  a  sequence  of  edge  flips.  Let  Do  and 
D\  be  the  restricted  Delaunay  triangulations  immediately  before  and  after  removing  a 
point  b.  We  argue  that  the  sequence  of  edge  flips  following  the  removal  successfully 
constructs  D\ . 

As  explained  in  Section  6,  only  the  polygon  bounded  by  the  link  of  b  requires  retrian¬ 
gulation.  The  algorithm  flips  one  diagonal  and  recurses  for  the  remaining  star  of  b  until 
only  three  triangles  remain.  It  thus  halts  after  a  number  of  edge  flips  that  is  less  than  the 
number  of  triangles  in  the  initial  star.  To  see  that  the  algorithm  is  correct,  we  observe 
that  each  flip  generates  a  triangle  that  is  guaranteed  to  belong  to  D\ .  The  membership  in 
D\  is  guaranteed  by  Function  IND,  which  checks  all  remaining  vertices  of  the  polygon 
and  not  just  one  as  for  the  flips  following  a  point  addition.  Edge  flips  that  cannot  be 
executed  because  one  of  the  endpoints  has  degree  3  will  again  not  be  attempted  because 
they  contradict  the  closed  ball  property  of  the  restricted  Voronoi  diagram. 


10.  Metamorphoses 

This  section  analyzes  the  point  configurations  generated  by  special  sampling.  Recall  that 
Hh  <  H  is  the  length  scale  threshold  that  triggers  the  start  and  end  of  special  sampling. 
In  the  forward  direction  we  start  with  a  two-sheeted  hyperboloid  that  enters  the  ball  with 
radius  Hh  around  its  center,  and  we  end  with  a  one- sheeted  hyperboloid  that  exits  the 
same  ball.  In  the  backwards  direction  the  events  are  the  same  in  reverse  order. 

Sizes  at  Transition.  Refer  to  the  double-cup  shown  in  Fig.  12.  The  t  4-  1  points  on  one 
sheet  form  a  regular  l-sided  cup.  The  i  vertices  of  the  base  lie  on  the  hot  circle  with 
radius  R0  =  Hyf  (  1  —  h2)/ 2,  which  lies  in  a  plane  at  distance  R\  =  H^/{  1  4-  h2)/ 2 
from  the  center.  Note  that  R%  +  R*  =  H 2.  Define  b  =  b,  and  c  =  bt+ 1,  with  indices 
modulo  i.  Independent  of  the  index  /,  the  lengths  of  the  edges  of  abc  are 

2 Rab  =  JWiiRi  -  Hh), 

I'Rbc  =  2/?o  sin  — . 

Any  isosceles  triangle  with  sides  of  length  E  and  height  L  has  circumradius  £2/2L. 

The  height  of  abc  is  La  bc  =  Rlb  —  R%c.  The  circumradius  is  therefore  AR^bl2Labc, 
which  is 

D  __  RARi-Hh) 

Kobe  —  ~~i .  - .  =. 

y/2Rl(Rl  -  Hh)  -  R2sin2(7T/£) 

Next  refer  to  the  cylinder- with-a- waist  shown  in  Fig.  12.  The  3 m  points  form  three 
parallel  regular  m-gons.  The  distance  between  two  contiguous  planes  is  Rq ,  and  the 
circumradii  of  the  three  m-gons  are  RuHh,R\.  Define  u  =  uif  v  =  m,*+ i,  w  =  u?if 
x  =  tUj.fi,  with  indices  modulo  m.  Independent  of  the  index  /,  the  lengths  of  the  edges 
uv  and  wx  are 

2RUV  =  2Hh  sin 

m 
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Fig.  18.  Portion  of  cylinder  in  Fig.  12  projected  onto  plane  parallel  to  m-gons. 


2  Rwx  =  2/?isin— . 

m 

To  compute  Rvw  ,RUVw>  Rvwx ,  we  consider  the  projection  of  the  middle  and  outer  m-gons 
onto  a  plane  parallel  to  the  two  m-gons,  as  shown  in  Fig.  18.  The  distance  between  the 
projections  of  w  and  (w  +  v)/2  is  Ri  —  Hh  cos(7r/m),  and  that  between  the  projections 
of  v  and  (w  +  x)/2  is  Ri  cos(;r/m)  -  Hh .  We  get  the  heights  LWyUV  and  Lv,wx  of  the 
two  triangles  by  taking  the  distances  to  three  dimensions,  which  means  squaring,  adding 
and  taking  square  roots.  The  length  of  an  edge  connecting  the  middle  m-gon  with 
one  of  the  two  outer  m-gons  is  the  root  of  R*v  4-  L2W  UV,  which  is 


=  yJ2R\(Ri  —  Hh  cos(7r/m)). 


We  compute  the  circumradii  of  the  two  isosceles  triangles  again  from  their  edges  and 
heights.  In  particular,  the  circumradius  of  uvw  is  4Rlw/2LWMVi  and  that  of  vwx  is 


Ri(R\  -  Hhcos(7t/m)) 

IH1  -  2R\Hh  cos(7r/m)  +  H2h2  cos2(7r/m) 
Ri(R\  —  Hhcos(7T/m)) 

I R]  ~~  2Ri  Hh  cos(jr/m)  +  R\  cos2(7r/m) 


Smooth  Transition.  We  derive  necessary  and  sufficient  conditions  for  h,  l,  m  that  guar¬ 
antee  a  smooth  transition  from  the  general  to  the  special  sampling  strategy.  By  this  we 
mean  that  the  configurations  at  the  start  of  a  metamorphosis  is  an  e  -sampling  and  satisfies 
both  Size  Bounds.  At  the  end  of  the  metamorphosis,  the  Size  Bounds  are  enforced  by 
eliminating  offending  edges  and  triangles  through  edge  contraction  and  vertex  insertion. 
The  result  is  a  triangulation  whose  vertex  set  is  an  ^-sampling  of  the  surface;  see  also 
the  remark  immediately  following  the  proof  of  the  Sampling  Lemma. 

The  length  scale  at  the  vertices  a,u,v  is  Hh,  and  that  at  b,  c,  w,  x  is  H.  The  Lower  and 
Upper  Size  Bounds  are  therefore  equivalent  to  Rbc  >  Ruvf  h ,  Rwx  ?  Rvw  >  (C/G)H 
and  Rabc,  Ruvw,  Rvwx  <  CQHh.  The  inequalities  for  Ruw,  Ruv,  Ruvw  are  redundant 
because  Rab  <  RVw,  Rwx  <  RUv/h,  RUVw  <  Rvwx  for  all  h  <  1.  In  addition  to  requiring 
that  the  triangles  abc,  uvw,  vwx  satisfy  the  Upper  Size  Bound,  it  is  convenient  to  require 
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also  that  their  radii  are  less  than  the  locally  allowed  minimum  edge  length.  This  extra 
requirement  implies  that  after  adding  points  on  and  inside  the  hot  sphere,  all  old  points 
inside  or  on  the  hot  sphere  are  too  close  to  at  least  one  new  point  and  thus  get  deleted. 
It  follows  that  all  remaining  old  vertices  lie  outside  the  hot  sphere.  We  thus  have  the 
following  two  conditions: 

(IV)  Rab/Hy  Rbc/H ,  Rwx/H  >  C/Q , 

(V)  Rabc/H,  Rvwx/H  <mm{Q,2/Q}Ch. 

Conditions  (I)-(V)  are  satisfied  for  e  =  0.279,  C  =  0.08,  Q  =  1.65,  h  =  0.98,  t  =  5, 
m  =  40.  We  summarize  the  results  assuming  this  assignment  of  constants. 

Transition  Lemma.  The  triangulation  at  the  start  of  a  metamorphosis  satisfies  the 
two  Size  Bounds  and  its  vertex  set  is  an  e -sampling  of  the  skin  surface . 

As  mentioned  earlier,  the  same  does  not  automatically  hold  for  the  end  configurations 
of  metamorphoses,  but  it  can  be  enforced  algorithmically.  The  purpose  of  bounding  the 
size  of  triangles  in  Condition  (V)  by  ( 2/Q)Ch  is  to  guarantee  that  the  algorithm  given 
in  Section  7  constructs  the  special  configurations  without  having  to  search  for  remaining 
old  vertices  inside  the  hot  sphere.  To  prove  this  algorithm  correct,  we  also  need  to  show 
that  these  configurations  are  part  of  the  restricted  Delaunay  triangulation,  which  follows 
from  the  Persistence  Lemmas  proved  below. 

Persistence  of  Triangulation.  We  show  that  the  special  configurations  exist  as  sub¬ 
complexes  of  the  restricted  Delaunay  triangulation  during  the  entire  time  interval  of  a 
metamorphosis.  Consider  the  simplices  in  the  Delaunay  complex  spanned  by  hot  vertices 
and  their  dual  Voronoi  polyhedra.  Fig.  1 9  sketches  both  for  the  double-cup  before  and  the 
cylinder- with-a-waist  after  the  middle  of  the  time  interval.  During  the  first  half  of  the  time 
interval,  the  hot  vertices  span  two  pyramids,  one  being  the  reflection  of  the  other  across  the 
symmetry  plane  of  the  hyperboloid.  The  points  are  in  degenerate  position,  which  implies 
that  the  Delaunay  complex  Ai  of  the  hot  points  contains  polyhedra  that  are  more  com- 


Fig.  19.  Solid  restricted  Delaunay  triangulations  and  dotted  Voronoi  polyhedra  of  hot  ball  configuration 
before  and  after  the  double-cone. 
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plicated  than  tetrahedra.  Specifically,  Ai  consists  of  the  two  pyramids  joined  by  an  edge 
connecting  their  apices  and  a  ring  of  five-sided  polyhedra  and  quadrangles  around  that 
edge.  As  usual,  D  denotes  the  restricted  Delaunay  triangulation  of  the  entire  vertex  set  V. 

Persistence  Lemma  A.  At  any  time  in  [-2 H2h2,  0),  the  intersection  of  Ai  and  D 
consists  of  the  two  rings  of  triangles  forming  the  double-cup  and  of  their  edges  and 
vertices. 

Proof  We  first  show  that  the  edges,  polygons,  and  polyhedra  in  Ai  that  do  not  belong 
to  the  double-cup  also  do  not  belong  to  D.  The  edges  connecting  the  two  cups  have 
dual  Voronoi  polygons  which  lie  in  the  symmetry  plane  separating  the  two  sheets  and 
therefore  cannot  intersect  the  hyperboloid.  To  see  that  they  do  not  intersect  any  other  part 
of  F ,  we  consider  the  sandwiching  spheres  defined  for  points  of  F  inside  the  hot  sphere. 
The  Voronoi  polygons  are  contained  in  the  union  of  balls  bounded  by  these  spheres,  else 
they  would  imply  an  empty  sphere  that  intersects  the  hyperboloid  in  a  patch  outside  the 
hot  sphere  that  is  large  enough  to  contradict  the  ^-sampling  property.  Detailed  compu¬ 
tations  of  a  lower  bound  for  the  size  of  such  an  implied  patch  are  omitted.  Since  D  is  a 
complex,  it  also  does  not  contain  the  Delaunay  polygons  and  polyhedra  incident  to  the 
excluded  edges.  The  base  polygons  of  the  two  pyramids  in  Ai  have  their  dual  Voronoi 
edges  on  the  symmetry  line  of  the  hyperboloid.  For  the  same  reason  as  above,  these 
edges  are  contained  in  the  union  of  balls  bounded  by  the  sandwiching  spheres  of  points 
of  F  inside  the  hot  sphere. 

We  second  show  that  the  triangles  abc  of  the  double-cup  belong  to  D .  At  time 
to  =  —2 H2h2  this  is  true  because  these  triangles  have  circumspheres  that  are  small 
enough  that  every  point  of  F  inside  these  spheres  would  belong  to  edges  that  violate  the 
Lower  Size  Bound.  At  times  t0  <  t  <  0  this  is  true  because  any  violation  is  prevented 
by  the  algorithm  before  it  occurs.  □ 

During  the  second  half  of  the  time  interval,  the  hot  vertices  form  three  convex  polygons 
in  three  parallel  planes.  The  middle  polygon  is  a  regular  m-g on  in  the  symmetry  plane 
of  the  hyperboloid,  and  the  other  two  are  reflections  of  each  other  across  that  plane  and 
are  inscribed  in  the  two  hot  circles.  The  Delaunay  complex  A2  of  the  hot  points  is  again 
degenerate,  consisting  of  the  above  mentioned  three  polygons,  which  form  the  top  and 
bottom  facets  of  two  drum-like  polyhedra.  The  two  drums  are  surrounded  by  a  ring  of 
four-sided  pyramids  alternating  with  tetrahedra. 

Persistence  Lemma  B.  At  any  time  in  (0,  2 H2h2],  the  intersection  of  A2  and  D  con¬ 
sists  of  the  two  rings  of  triangles  forming  the  cylinder-with-a-waist  and  of  their  edges 
and  vertices. 

Proof  The  edges,  polygons,  and  polyhedra  in  A2  that  do  not  belong  to  the  cylinder- 
with-a-waist  have  their  dual  Voronoi  polygons,  edges,  and  vertices  either  in  the  symmetry 
plane  or  the  symmetry  axis  of  the  hyperboloid.  For  the  reason  mentioned  in  the  proof 
of  the  Persistence  Lemma  A,  these  polygons,  edges,  and  vertices  are  contained  in  the 
union  of  balls  bounded  by  sandwiching  spheres  of  points  of  F  inside  the  hot  sphere.  The 
corresponding  edges,  polygons,  and  polyhedra  of  A2  thus  do  not  belong  to  D. 
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The  remainder  of  the  proof  establishes  that  the  triangles  uvw  and  vwx  belong  to 
D .  Immediately  after  time  /  =  0  this  is  true  because  the  triangles  in  the  double-cup 
belonged  to  D  immediately  before  time  t  =  0.  At  times  0  <  t  <  t\  this  is  true  because 
any  violation  is  prevented  by  the  algorithm  before  it  occurs.  □ 

The  two  Persistence  Lemmas  also  hold  for  the  reverse  metamorphosis,  which  changes 
a  one-sheeted  into  a  two-sheeted  hyperboloid.  To  see  this,  run  time  backwards  and  ex¬ 
change  the  arguments  that  establish  that  the  two  special  configurations  are  subcomplexes 
of  D  when  they  are  first  created.  These  arguments  are  contained  in  the  respective  last 
paragraphs  of  the  two  proofs. 

Summary.  The  two  Persistence  Lemmas  establish  that  the  closed  ball  property  of  the 
restricted  Voronoi  diagram  is  maintained  even  inside  the  hot  spheres  that  guide  the 
algorithm  through  the  various  metamorphoses. 

Special  Homeomorphism  Theorem.  The  restricted  Delaunay  triangulation  of  the 
points  chosen  by  special  sampling  triangulates  the  skin  surface  inside  each  hot  sphere. 

Together  with  the  General  Homeomorphism  Theorem  this  implies  that  we  have  a 
triangulation  of  the  skin  surface  at  all  times. 


11.  Discussion 

This  paper  describes  a  dynamic  algorithm  for  maintaining  the  triangulation  of  a  deform¬ 
ing  skin  surface  by  adapting  it  to  changing  shape,  curvature,  and  topology. 

Abstract  Interface.  The  algorithm  uses  detailed  knowledge  about  the  skin  surface  to 
avoid  pitfalls,  such  as  insufficient  quality  of  approximation,  small  angles,  and  wrong 
connections.  In  an  effort  to  understand  the  extent  to  which  the  algorithm  can  be  general¬ 
ized,  we  may  ask  how  much  knowledge  about  the  surface  the  algorithm  really  needs.  Can 
we  list  axioms  for  a  deforming  surface  that  imply  the  applicability  of  the  algorithm?  To 
make  this  a  worthwhile  exercise,  one  would  of  course  hope  that  the  class  of  surfaces  and 
deformations  defined  by  the  axioms  is  significantly  larger  than  the  class  of  skin  surfaces 
and  the  growth  model. 

A  key  idea  of  the  algorithm  is  to  keep  the  density  of  vertices  roughly  proportional  to 
one  over  the  maximum  curvature.  This  is  only  possible  if  that  measure  of  curvature 
satisfies  a  one-sided  Lipschitz  condition,  like  that  stated  in  the  Curvature  Variation 
Lemma.  It  is  conceivable  that  the  sampling  density  can  be  based  on  other  expressions  of 
curvature  that  satisfy  a  one-sided  Lipschitz  condition.  The  local  feature  size  proposed  in 
[3]  is  a  candidate  for  such  an  expression,  but  it  is  usually  not  easy  to  compute.  Another 
important  cornerstone  of  the  algorithm  is  the  predictability  of  metamorphoses.  Without 
knowing  when  and  where  a  metamorphosis  will  happen,  we  would  have  to  resort  to 
relatively  expensive  collision  detection  and  resolution  algorithms. 

Implementation.  Where  do  we  go  from  here?  The  first  job  on  the  agenda  is  the  rigorous 
implementation  of  the  dynamic  skin  triangulation  algorithm.  The  first  author  of  this  paper 
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has  already  taken  steps  in  that  direction,  partially  reusing  prior  software  on  alpha  shapes 
[9]  and  on  computing  Betti  numbers  [6].  It  will  be  interesting  to  study  the  algorithm 
experimentally  and  measure  the  influence  of  design  decisions  on  its  performance.  For 
example,  the  bounds  on  the  constants  C,  Q  controlling  the  curvature  adaptation  algorithm 
derived  in  this  paper  are  all  conservative.  Perhaps  it  is  possible  to  relax  the  requirements 
a  fair  amount  without  compromising  the  correctness  of  the  algorithm. 

It  would  be  interesting  to  modify  the  algorithm  to  sample  points  with  local  density 
roughly  proportional  to  *JTJk  rather  than  to  1/k.  A  density  like  that  would  lead  to  a 
uniformly  approximating  triangulation.  In  other  words,  the  distance  between  the  trian¬ 
gulation  and  the  skin  surface  would  be  bounded  by  a  constant  independent  of  curvature. 
As  another  bonus,  the  vertices  would  be  distributed  without  accumulation  points  at  sin¬ 
gularities  with  infinite  curvature.  As  a  corresponding  drawback,  the  triangulation  would 
suffer  from  violations  of  the  closed  ball  property,  which  is  crucial  to  the  algorithm  as 
described  in  this  paper. 

Deformation.  Recall  that  the  skin  surface  is  defined  by  a  finite  set  of  spheres  in  R3.  It 
is  deformed  by  manipulating  this  set.  The  simplest  type  of  deformation  is  plain  growth, 
which  is  generated  by  increasing  all  radii  continuously  and  simultaneously  in  a  way  that 
preserves  the  mixed  complex.  The  next  more  complicated  type  of  deformation  is  the  one 
described  in  [4].  It  can  be  used  to  morph  one  skin  surface  into  any  other  skin  surface 
by  a  combination  of  moving  and  growing/ shrinking  of  spheres.  This  is  done  in  a  way 
that  preserves  the  combinatorics  but  not  the  geometric  realization  of  the  mixed  complex. 
Metamorphoses  in  the  morphing  model  are  only  slightly  more  difficult  to  predict  than 
in  the  growth  model.  The  maintenance  of  the  skin  triangulation  is  however  made  more 
complicated  by  the  occurrence  of  new  types  of  degenerate  metamorphoses,  such  as  the 
ones  that  progress  to  the  critical  point  and  are  then  reversed.  In  any  physical  simulation, 
where  the  deformation  depends  on  forces  unrelated  to  the  combinatorics  of  the  mixed 
complex,  we  will  also  have  to  maintain  that  complex  dynamically.  The  computational 
overhead  is  negligible  if  we  use  a  dynamic  algorithm  for  maintaining  the  Delaunay 
complex  [11].  Note  that  a  single  time  step  in  the  simulation  may  jump  over  any  number 
of  metamorphoses  and  other  changes.  A  reasonable  approach  to  bridging  the  gaps  in 
time  is  to  connect  any  two  contiguous  time  slices  with  a  deformation  in  the  morphing 
model  mentioned  above. 

An  important  question  in  this  context  is  the  inverse  problem.  How  do  we  construct 
the  skin  surface  that  best  approximates  a  given  shape,  and  how  do  we  maintain  such 
an  approximation?  We  need  modeling  operations  that  allow  us  to  increase  or  decrease 
the  complexity  of  the  surface  locally.  The  former  can  be  done  by  doubling  spheres  and 
continuously  moving  them  apart.  The  latter  can  be  done  by  the  inverse  operation,  which 
moves  spheres  together  and  merges  them  into  one. 
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Appendix 

Figure  20  gives  an  overview  of  the  proof  structure  used  in  this  paper.  Nodes  represent 
properties,  lemmas,  invariants,  and  conditions.  Arcs  represent  logical  dependencies  di¬ 
rected  from  top  to  bottom.  All  abbreviations  are  explained  below.  Table  2  provides  a  list 
of  notation  used  in  this  paper. 


0) 

Condition  on  e 

Section  4 

(D) 

Condition  on  C  and  Q 

Section  8 

ail) 

Condition  on  C  and  Q 

Section  8 

(IV) 

Condition  onh,Z,m 

Section  10 

(V) 

Condition  on  h,  l,  m 

Section  10 

[L] 

Lower  Size  Bound 

Section  6 

[U] 

Upper  Size  Bound 

Section  6 

CcL 

Circumcenter  Lemma 

Section  8 

CSL 

Curvature  Sandwich  Lemma 

Section  3 

CVL 

Curvature  Variation  Lemma 

Section  3 

ENL 

Edge  Normal  Lemma 

Section  4 

GHT 

General  Homeomorphism  Theorem 

Section  4 

HSL 

Hot  Spot  Lemma 

Section  7 

IcL 

Iso-curvature  Lemma 

Section  3 

LDC 

Long  Distance  Claim 

Section  4 

MAL 

Minimum  Angle  Lemma 

Section  6 

NSL 

No-Short-Edge  Lemma 

Section  8 
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Fig.  20.  Logical  structure  of  the  proof  that  the  algorithm  in  this  paper  constructs  a  triangulation  of  the  skin 
surface. 
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Table  2.  Notation  for  important  geometric  concepts,  functions,  variables, 
and  constants. 


Xl.X2.X3 

Coordinates 

7 Tq.  R3  — >•  R 

Weighted  (square)  distance  function 

a  =  (a.  A) 

Zero-set  of  7r£,  sphere  with  center  a  and  radius  A 

aa  +  fib 

Zero- set  of  an^  + 

VI  =  (a,  A/V 2) 

Zero-set  of  —  7t^(a)fl 

F  —  skin  A 

Skin  surface,  env  V conv  A 

body  A 

Body  bounded  by  skin  A 

vx 

Intersection  of  Voronoi  polyhedra 

Sx 

Delaunay  simplex 

k 

Dimension  of  Delaunay  simplex 

fix 

Mixed  cell,  ( vx  +  ^a,)/2 

zx 

Center,  aff  vx  D  aff  8x 

k:  F  ->  R 

Maximum  curvature  function 

p:  F  ->  R 

Length  scale  function,  \/k 

nx » nxyz 

Normal  vector  at  a:,  of  xyz 

Tangent  vector,  (y  -  x)/\\y  -  a|| 

Sandwiching  spheres  at  jc  e  F 

Rabc 

Radius  or  size  of  edge,  triangle 

Qab 

Target  size,  max{g(a),  g(b)} 

Qabc 

Target  size,  min{p(a),  g{b).  p(c)} 

£ 

Constant  controlling  sampling  density 

Q 

Constant  controlling  triangle  quality 

c 

Constant  controlling  approximation 

M:  R3  R 

Trajectory  of  skin  surface 

<Pij 

Diffeomorphism  between  skins 

V 

^-Sampling  of  F 

D 

Restricted  Delaunay  triangulation 

A 

Delaunay  complex 

H 

Radius  of  hot  ball 

Hot  portion  of  space 

Fh 

Hot  portion  of  skin  surface,  F  fi  R^ 

h 

Constant  triggering  a  metamorphosis 

l.m 

Constant  numbers  of  hot  vertices 

to.h 

Start /end  time  of  a  metamorphosis 

Ro.  Ri 

Start/end  radius  of  hot  circles 

NVL 

Normal  Variation  Lemma 

Section  3 

PLA 

Persistence  Lemma  A 

Section  10 

PLB 

Persistence  Lemma  B 

Section  10 

SDC 

Short  Distance  Claim 

Section  4 

SgL 

Sampling  Lemma 

Section  8 

SHT 

Special  Homeomorphism  Theorem 

Section  10 

SwP 

Sandwich  Property 

Section  2 

TNL 

Triangle  Normal  Lemma 

Section  4 

TsL 

Transition  Lemma 

Section  10 

VEL 

Voronoi  Edge  Lemma 

Section  4 

VGL 

Voronoi  Polygon  Lemma 

Section  4 

VHL 

Voronoi  Polyhedron  Lemma 

Section  4 
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